corrected length of lines, ran with opty-dev, only two obj functions …

…considered for speed

examples-gallery/plot_sliding_block.py

@@ -4,7 +4,7 @@
 =========================
 A block, modeled as a particle is sliding on a road to cross a
 hill. The block is subject to gravity and speed dependent
-friction. 
+friction.
 Gravity points in the negative Y direction.
 A force tangential to the road is applied to the block.
 Three objective functions to me minimized may be selected:
@@ -36,10 +36,9 @@
 import sympy as sm
 from opty.direct_collocation import Problem
 import matplotlib.pyplot as plt
+import matplotlib as mp
 from matplotlib.animation import FuncAnimation
 from copy import deepcopy
-
-
 # %%
 # define the road
 def strasse(x, a, b):
@@ -74,10 +73,10 @@ def strasse(x, a, b):
 
 q_ind = [x]
 u_ind = [ux]
- 
+
 KM = me.KanesMethod(N, q_ind=q_ind, u_ind=u_ind, kd_eqs=kd)
 (fr, frstar) = KM.kanes_equations(BODY, FL)
-EOM = kd.col_join(fr + frstar) 
+EOM = kd.col_join(fr + frstar)
 EOM.simplify()
 sm.pprint(EOM)
 
@@ -87,7 +86,7 @@ def strasse(x, a, b):
 speicher = deepcopy((x, ux, F))
 
 # which objective functions are to be used
-selektor = (2, 1)
+selektor = (1, 0)
 
 # store results of the optimizations for later plotting
 solution_list = [0., 0., 0.]
@@ -125,11 +124,11 @@ def strasse(x, a, b):
         interval_value = h
 
     if selektion == 1:
-        def obj(free): 
-            Fx = free[laenge * num_nodes: (laenge + 1) * num_nodes] 
+        def obj(free):
+            Fx = free[laenge * num_nodes: (laenge + 1) * num_nodes]
             return free[-1] * np.sum(Fx**2) + free[-1] * verstaerkung
 
-        def obj_grad(free):  
+        def obj_grad(free):
             grad = np.zeros_like(free)
             l1 = laenge * num_nodes
             l2 = (laenge + 1) * num_nodes
@@ -145,10 +144,10 @@ def obj_grad(free):
             grad = np.zeros_like(free)
             grad[-1] = 1.
             return grad
-    
+
     elif selektion == 2:
-        def obj(free): 
-            Fx = free[laenge * num_nodes: (laenge + 1) * num_nodes] 
+        def obj(free):
+            Fx = free[laenge * num_nodes: (laenge + 1) * num_nodes]
             return interval_value * np.sum(Fx**2)
 
         def obj_grad(free):
@@ -164,51 +163,53 @@ def obj_grad(free):
 
 # pick the integration method. backward euler and midpoint are the choices.
 # backward euler is the default.
-    methode = 'backward euler' 
+    methode = 'backward euler'
 
-    if selektion in (0, 1): 
-        initial_guess = np.array(list(np.ones((len(state_symbols) 
+    if selektion in (0, 1):
+        initial_guess = np.array(list(np.ones((len(state_symbols)
             + len(specified_symbols)) * num_nodes) * 0.01) + [0.02])
     else:
-        initial_guess = np.ones((len(state_symbols) + 
+        initial_guess = np.ones((len(state_symbols) +
             len(specified_symbols)) * num_nodes) * 0.01
 
     initial_state_constraints = {x: 0., ux: 0.}
 
     final_state_constraints = {x: 10., ux: 0.}
-
-    instance_constraints = ( 
-        tuple(xi.subs({t: t0}) - xi_val for xi, xi_val in initial_state_constraints.items()) + 
-        tuple(xi.subs({t: tf}) - xi_val for xi, xi_val in final_state_constraints.items())
+
+    instance_constraints = (
+        tuple(xi.subs({t: t0}) - xi_val for xi, xi_val
+            in initial_state_constraints.items()) +
+        tuple(xi.subs({t: tf}) - xi_val for xi, xi_val
+            in final_state_constraints.items())
     )
 
 # forcing h > 0 avoids negative h as 'solutions'.
     if selektion in (0, 1):
-        bounds = {F: (-15., 15.), x: (initial_state_constraints[x], 
+        bounds = {F: (-15., 15.), x: (initial_state_constraints[x],
         final_state_constraints[x]), ux: (0., 1000.), h:(1.e-5, 1.)}
     else:
-        bounds = {F: (-15., 15.), x: (initial_state_constraints[x], 
+        bounds = {F: (-15., 15.), x: (initial_state_constraints[x],
         final_state_constraints[x]), ux: (0., 100)}
 
-    prob = Problem(obj, 
-        obj_grad, 
-        EOM, 
-        state_symbols, 
-        num_nodes, 
+    prob = Problem(obj,
+        obj_grad,
+        EOM,
+        state_symbols,
+        num_nodes,
         interval_value,
         known_parameter_map=par_map,
         instance_constraints=instance_constraints,
         bounds=bounds,
         integration_method=methode)
 
 # set max number of iterations. Default is 3000.
-    prob.add_option('max_iter', 3000)          
+    prob.add_option('max_iter', 3000)
 
     solution, info = prob.solve(initial_guess)
     solution_list[selektion] = solution
     info_list[selektion] = info
     prob_list[selektion] = prob
-# %% 
+# %%
 # animate the solutions and plot the results
 def drucken(selektion, fig, ax, full=True):
     solution = solution_list[selektion]
@@ -223,15 +224,15 @@ def drucken(selektion, fig, ax, full=True):
     interval_value = duration / (num_nodes - 1)
 
     strasse1 = strasse(x, a, b)
-    strasse_lam = sm.lambdify((x, a, b), strasse1, cse=True)    
+    strasse_lam = sm.lambdify((x, a, b), strasse1, cse=True)
 
     P0_x = solution[:num_nodes]
     P0_y = strasse_lam(P0_x, par_map[a], par_map[b])
 
 # find the force vector applied to the block
     alpha = sm.atan(strasse(x, a, b).diff(x))
     Pfeil = [F*sm.cos(alpha),  F*sm.sin(alpha)]
-    Pfeil_lam = sm.lambdify((x, F, a, b), Pfeil, cse=True)     
+    Pfeil_lam = sm.lambdify((x, F, a, b), Pfeil, cse=True)
 
     l1 = laenge * num_nodes
     l2 = (laenge + 1) * num_nodes
@@ -246,10 +247,10 @@ def drucken(selektion, fig, ax, full=True):
 
     if full == True:
         print('message from optimizer:', info['status_msg'])
-        print(f'objective value {obj(solution):,.1f}')
         if selektion in (0, 1):
             print(f'optimal h value is: {solution[-1]:.3f}')
         prob.plot_objective_value()
+        prob.plot_constraint_violations(solution)
         prob.plot_trajectories(solution)
 
     def initialize_plot():
@@ -262,40 +263,44 @@ def initialize_plot():
         if selektion == 0:
             msg = f'speed was optimized'
         elif selektion == 1:
-            msg = f'speed and energy optimized, weight of speed = {verstaerkung:.1e}'
+            msg = (f'speed and energy optimized, weight of speed = ' +
+                f'{verstaerkung:.1e}')
         else:
             msg = f'energy optimized'
 
         ax.grid()
         strasse_x = np.linspace(xmin, xmax, 100)
         ax.plot(strasse_x, strasse_lam(strasse_x, par_map[a], par_map[b]),
             color='black', linestyle='-', linewidth=1)
-        ax.axvline(initial_state_constraints[x], color='r', linestyle='--', linewidth=1)
-        ax.axvline(final_state_constraints[x], color='green', linestyle='--', linewidth=1)
+        ax.axvline(initial_state_constraints[x], color='r',
+            linestyle='--', linewidth=1)
+        ax.axvline(final_state_constraints[x], color='green',
+            linestyle='--', linewidth=1)
 
 # Initialize the block and the arrow
         line1, = ax.plot([], [], color='blue', marker='o', markersize=12)
-        pfeil   = ax.quiver([], [], [], [], color='green', scale=35, width=0.004)
+        pfeil   = ax.quiver([], [], [], [], color='green', scale=35,
+            width=0.004)
         return line1, pfeil, msg
 
     line1, pfeil, msg = initialize_plot()
 
 # Function to update the plot for each animation frame
     def update(frame):
-        message = (f'running time {times[frame]:.2f} sec \n' +  
-            f'the red line is the initial position, the green line is the final position \n' + 
+        message = (f'running time {times[frame]:.2f} sec \n' +
+            f'the red line is the initial position, the green line is ' +
+            f'the final position \n' +
             f'the green arrow is the force acting on the block \n' +
             f'{msg}' )
         ax.set_title(message, fontsize=12)
-    
+
         line1.set_data([P0_x[frame]], [P0_y[frame]])
         pfeil.set_offsets([P0_x[frame], P0_y[frame]])
         pfeil.set_UVC(Pfeil_x[frame], Pfeil_y[frame])
         return line1, pfeil
 
-    animation = FuncAnimation(fig, update, frames=range(len(P0_x)), 
+    animation = FuncAnimation(fig, update, frames=range(len(P0_x)),
         interval=1000*interval_value, blit=True)
-
     return animation, update
 
 ## %%
@@ -322,7 +327,7 @@ def update(frame):
 # A frame from the animation.
 fig, ax = plt.subplots(figsize=(8, 8))
 _, update = drucken(selektor[0], fig, ax, full=False)
-# sphinx_gallery_thumbnail_number = 5
+number = len(selektor) * 4 + 1
+# sphinx_gallery_thumbnail_number = number
 update(100 )
-
 plt.show()