A question: Form of the EOMs for opty.

[Peter230655]

I am fairly familiar with sympy.physics.mechanics, but totally new to opty. When I create the symbolic EOMs with symy mechanics, I always do it löike this:

KM = me.KanesMethod(N, q_ind=q_ind, u_ind=u_ind, kd_eqs=kd) (fr, frstar) = KM.kanes_equations(BODY, FL) MM = KM.mass_matrix_full force = KM.forcing_full

and I solve for $\dot y = MM^{-1} \cdot force$ only numerically during the integration. If I understood the documentation correctly, here I should do like this: RHS = KM.rhs() #This symbolically forms $MM^{-1} \cdot force$, I believe. I normally do not use it, as the result may be a 'monster'

EOM = $\dot y - RHS$

Is my thinking correct? Thanks!

[moorepants]

No, that is a waste of symbolic computation. opty only needs the equations in the form fr + frstar for an unconstrained system. For a constrained system you also need to include the holonomic constraints. Solving the linear system is not necessary, but if you have equations in that form you can do the ydot - rhs as you mention.

[tjstienstra]

Here is also an awesome video, which I remember watching to also learn some of the basics about direct collocation and its usage for solving trajectory optimization problems. I should note that he also discusses some stuff which is beyond opty (opty is relatively simple), but the video is in general not too difficult to follow. I remember also reading some other stuff, but I haven't been able to immediately find that

[Peter230655]

Thanks! so I would do like this (I will start with non constraint systems - better for a beginner, I think):

EOM = sm.Matrix.vstack([fr, frstar]

If correct could I do the same, when I use solve_ivp?

@Timo: I will watch this video!

[tjstienstra]

Yes, definitely start with unconstrained systems.

Remember that fr and frstar are the generalized active forces and generalized inertial forces, so you need to sum them to get the equations of motion: eoms = fr + frstar

As the video I think also explains. Direct collocation treats the eoms as constraints, so they should be zero, while solve_ivp is an integrator requiring the gradient of the function you'd like to integrate, i.e. speed and acceleration.

[Peter230655]

In the middle of the vidoe, really good!

Thanks! I never gave much thought to fr, frstart. Need to look at Jason's lecture again!

This video is really good!!

[Peter230655]

Sorry that I still have to come back to the issue of the OEMs for opty. Eq(186) of Jason's lecture sates: $\bar F_r + \bar F^{\star}_r = f_d(\dot{\bar u}, \bar u, \bar q, t) = 0$ So, if $\bar{kd}$ are my kinematic equations, I would have to write: EOM = sm.Matrix.vstack([ $\bar{kd}, \bar F_r + \bar F^{\star}_r $]) Would this be right? Would this be what Timo worte in line 25 of cell 2 in his program:

eoms = kdes.col_join(Fr + Frs)

Thanks for any help and clarification!

[moorepants]

Yes, you need to stack the kinematical and dynamical equations. You can see how that looks here: https://github.com/csu-hmc/opty/blob/master/examples/pendulum_swing_up.py#L35 for a SDOF pendulum.

[Peter230655]

Thanks! This is the example I am looking at right now, Timo's example is more difficult. Now clear oth this point.

[Peter230655]

As reference, I use this: https://opty.readthedocs.io/_/downloads/en/latest/pdf/ I am still playing with the inverted pendulum. This program creates a number of plots, tow of which I do not understand: 1. Amongst the plots comstraint violations, the third one is Instance. I do not understand what this means. 2. There is a plot objective value. Initially, I thought this was the value of the cost function - but this cannot be the case, as it 'goes up and down', while the cost function $T(t)^2 \geq 0.$ Thanks for the help!!

[tjstienstra]

Amongst the plots comstraint violations, the third one is Instance. I do not understand what this means.

The plot constraints gives a measure of how much constraints are violated in the final solution. It shows the violation of the equations of motion. And the violation of the instance constraints. The instance constraints are constraints on what certain states should have at certain moments in time. For the inverted pendulum we state that the angle of the pendulum at t=0 should be 0 degrees, and at t=10.0 it should be π. Likewise the angular velocity should be 0 at both the start and the end. In this plot we can see that all are met (violation in the order of 1e-32 for both angular velocities).

There is a plot objective value. Initially, I thought this was the value of the cost function - but this cannot be the case, as it 'goes up and down', while the cost function

The objective value is plotted against the iterations of the interior point solver. This says something about the convergence of the problem.

[Peter230655]

Thanks! I think I understood. In the example program of the inverted pendulum (I want to understand this before I tackle Timo's example) there is this:

_num_nodes = 500

I, m, g, d, t = sym.symbols('I, m, g, d, t') theta, omega, T = sym.symbols('theta, omega, T', cls=sym.Function)

state_symbols = (theta(t), omega(t)) constant_symbols = (I, m, g, d) specifiedsymbols = (T(t),)

Then, further down is this: _solution, info = prob.solve(initialguess)

solution is an np.array of shape (1500, ). I assume, the first 500 entries are the values of theta at the node points, the next 500 entries are the values of omega at the node points, the final 500 entries are the values of T at the node points. Is this correct?

If correct, would solution sort like this: solution[: num_nodes] = values of state_symbols[0] solution[num_nodes: 2*num_nodes = values of state_symbols[1]

and after all the state_symbols are taken care of, the entries of specified_symbols would follow? Sorry about all these questions - and thans for the help!

[tjstienstra]

I assume, the first 500 entries are the values of theta at the node points, the next 500 entries are the values of omega at the node points, the final 500 entries are the values of T at the node points. Is this correct?

Yes, that is correct. The shape of the free array is shape(n*N + q*N + r, 1), where

• N : number of collocation nodes
• n : number of unknown state trajectories
• q : number of unknown input trajectories
• r : number of unknown parameters

So the first n N are states, the next q N are input trajectories, e.g. forces and the last r are unknown parameters, e.g. if you would like to optimize something for a constant like the mass.

If correct, would solution sort like this: solution[: num_nodes] = values of state_symbols[0] solution[num_nodes: 2*num_nodes = values of state_symbols[1]

Yes this is entirely correct.

and after all the state_symbols are taken care of, the entries of specified_symbols would follow?

Refer to the explanation above. Pretty sure that explanation is also somewhere in the documentation, I just copied parts of it from random locations in opty's source code.

[Peter230655]

Thanks! Source code would be like Chinese to me - and I would not be able to find in in the first place! :-)

I guess your confirmation above also explains why

_def obj(free): T = free[2 num_nodes:] return interval_value np.sum(T**2)

def obj_grad(free): grad = np.zeros_like(free) grad[2 num_nodes:] = 2.0 interval_value free[2 numnodes:] return grad

are put this way: The cost function and its gradient are only aspplied to T ?

[tjstienstra]

are put this way: The cost function and its gradient are only aspplied to T ?

Yes, because we want to minimize the torque. Last week I implemented a function in opty that actually generates these two functions based on an expression like sm.Integral(F(t)**2, t).

A bit simpler and primitive version was included in the file I sent you earlier.

[Peter230655]

Oh, so you added the featurew of creating symbolic cost functionals?! As I understood the latest documentation Jason wrote it this was not available. So now it is? I will study the file you sent only after I understand the inverted pendulum - it is simpler! :-)