Generate trapezoidal collocation constraints

[DeepaMahm]

Hi @MatthewPeterKelly ,

Excuse me fro the naive question. I'm new here
I'd like to know if the user should always define the constraints (eg. non-linear equality constraints)
for solving the trajectory optimization problem.
For example, The system dynamics that I want to provide as a set of constraints is of the form
dS/dt = -AS. (this represents the dynamics of a graph network)

S is a vector and A is a matrix.

I would like to know if it is possible to automatically generate trapezoidal collocation constraints using OptimTraj. I'd like to supply these constraints to fmincon

The complete description of the problem can be found here

[MatthewPeterKelly]

Hi @DeepaMahm -- short answer: Yes.

OptimTraj has trapezoidal collocation as one of the built-in methods. When you pass in the options struct, set the .method field to trapezoid. You need to tell the optimization about your system dynamics dS/dt = -A*S, which you do by creating a function that performs that calculation and then setting that function to the problem.func.dynamics field of the problem setup.

Good luck!

[DeepaMahm]

Hi @MatthewPeterKelly Thanks a lot for the response. I'm referring to the user guide for setting p the problem.
I'm trying to set up the field in problem.guess. The control parameters in my system are not a function of time and I am not sure how the following has to be set up.

problem.guess.control = ug size(ug) = [Nu, Ng]
In brief, I've the following equations
Equation1 : (represents the exact dynamics of a flow in a graph network)

dS/dt = - M^T D M S + other effects
Equation 2 : (represents the approximate version of equation 1)
dS_approx/dt = - M^T D_hat M S

In my objective function, I try to find a D_hat the minimizes the difference between the time course data predicted by equation 1 and 2.

So, D_hat is not essentially a function of time. (D_hat is a diagonal matrix and the diagonal entries are the edge weights of a graph)

For a detailed explanation of the above system, please check here

In such a case, I would like to know how
problem.guess.control = ug size(ug) = [Nu, Ng] has to be defined.

Second doubt:
I'm referring to line 113 of file trapezoid.m in OptimTraj
[defects, defectsGrad] = computeDefects(dt,x,f,dtGrad,xGrad,fGrad)
From the description given for this function,

% INPUTS:
% dt = time step (scalar)
% x = [nState, nTime] = state at each grid-point along the trajectory
% f = [nState, nTime] = dynamics of the state along the trajectory

If my understanding is right, x and f should be returned from equation 1
dS/dt = - M^T D M S + other effects
that represents the exact dynamics of my system.

Could you please clarify if x = S that is obtained after solving the ode in equation (1) and
f = - M^T D M S + other effects

computed at stored at each time step.

Also, I could understand that
nonlcon = % should call a function that returns [c ceq], in OptimTraj collectConstraints returns [c ceq]
for passing it to fmincon for parameter estimation
Dhat= fmincon(@objfun,Dhat0,A,b,Aeq,beq,lb,ub,nonlcon)
But , I'm not sure which funuction passes the input arguments ( dt, x and f) for computeDefects(dt,x,f,dtGrad,xGrad,fGrad) in OptimTraj

Could you please clarify?
I'm sorry for the really long post
Thanks a lot

[DeepaMahm]

Hello @MatthewPeterKelly
I'm really sorry for posting again, but I couldn't really understand how the input arguments are passed to function,
function [defects, defectsGrad] = computeDefects(dt,x,f,dtGrad,xGrad,fGrad) in line 113 of file trapezoid.m

Could you please explain?
Thanks a lot.

[MatthewPeterKelly]

Here are some quick answers that should help you make progress in the right direction:

1. What is going on with computeDefects? It is complicated because both trapezoid and hermiteSimpson are set up to use the same direct collocation code in the background. Read a bit more about how anonymous (lambda) functions work in Matlab (the @ symbol). You'll eventually find out that the computeDefects() function is actually called on line 291 of directCollocation.m. You should need to worry about these details though for your problem.

2. Without digging into the details, it sounds like you have a parameter optimization problem. In practice, this means that your "control" inputs happen to be constant. For example, suppose that your dynamics are x_dot = k * A * x, where k is an unknown scalar constant and A is a known constant matrix. You can pose this as a trajectory optimization problem by creating a new state for k that has a zero derivative no boundary constraints. I forget if my framework allows you to solve problems without a traditional control, but if not then you could fake it by making a dummy control that has a quadratic cost to push it to zero.

[MatthewPeterKelly]

Quick follow-up on 2: I have not tried this with OptimTraj, and may be forgetting some detail. The software GPOPS-II has good support for this type of optimization, and they clearly explain in their user guide (or at least used to) how it was implemented internally. I suggest that you look up some of their papers or at least read their user guide if you want to go down this route.

[DeepaMahm]

Hi @MatthewPeterKelly
Thanks for letting me know about GPOPS-II. However, I couldn't find a free license and I'm continuing to implement using OptimTraj.

A simplified explanation of the problem that I am trying to solve.

---

then you could fake it by making a dummy control that has a quadratic cost to push it to zero.

I've done this in the following function

function f = objfun(Dhat)

%% Integrator settings
phi0    = [5; 0; 0; 0; 0; 0; 0; 0; 0; 0];
tspan   = 0:dt:0.5;
options = odeset('abstol', 1e-10, 'reltol', 1e-9);

%% generate exact solution
    [t, phi]  = ode15s(@(t,phi) experimental(t,phi), tspan , phi0 ,options);


%% generate approximate solution

    [t, phi_tilde]  = ode15s(@(t,phi_tilde) predicted(t,phi_tilde, Dhat), tspan , phi0 ,options);


%% objective function for fminunc/fmincon
    f = (phi - phi_tilde).*(phi - phi_tilde);
    f = sum(f, 'all')
    
%% objective function for lsqnonlin
%     f  = phi - phi_tilde
%     f = sum(f, 'all')

end

I've a doubt in how

xLow = phi_tilde(:,idxLow);
xUpp = phi_tilde(:,idxUpp);

fLow = f(:,idxLow);
fUpp = f(:,idxUpp);

are computed in computeDefects. Could you please let me know if the
x = [nState, nTime] = state at each grid-point along the trajectory
f = [nState, nTime] = dynamics of the state along the trajectory

are pre-computed by solving the dynamics with control parameters (equation 2 in my system) using ode solver.
Something like below,

[t, x]  = ode15s(@(t,phi_tilde) predicted(t,phi_tilde, Dhat), tspan , phi0 ,options);
f =cellfun(@predicted,num2cell(t),num2cell(phi_tilde,2),num2cell(Dhat_mat,2),'UniformOutput',false);
 

Also, I would like to know if I can contribute this to OptiTraj demo examples if I can successfully set up and solve.

[DeepaMahm]

Hi @MatthewPeterKelly Sorry for bothering you. I tried to implement the defects in my problem, adapting the function available in OptimTraj. I could get the optimal solution without constraints, but with constraints, there is a problem. I am also attaching the code here for your kind reference.

[MatthewPeterKelly]

I quickly skimmed through the math of your problem statement. I see two issues with your implementation: 1) You should not be discretizing the problem yourself - that it the job of OptimTraj. Let it do the work for you. If you want to have the output of Eqn1 and Eqn2 match, then make each of them into a state variable, let OptimTraj do the dynamics, and then put a cost between them. As a general rule, if you are ever trying to guess how the discretization works in order to write your own objective or constraint, then you're setting something up wrong. 2) Calling any of Matlab's ode*() functions inside of a function handle that is passed into OptimTraj will cause poor results. This is because all of the ode*() functions are accurate but not consistent: they perform a different sequence of operations on each call. This causes discontinuities in gradients between iterations and dramatically slows down FMINCON. You can get rid of this "simulation" step by using the integration that is built into OptimTraj - it will make your code easier to write, the solution will be more accurate, and it will run much faster. Good luck! Matt

…

On Tue, Mar 31, 2020 at 11:38 AM DeepaMahm ***@***.***> wrote: Hi @MatthewPeterKelly <https://github.com/MatthewPeterKelly> Sorry for bothering you. I tried to implement the defects in my problem, adapting the function available in OptimTraj. I could get the optimal solution without constraints, but with constraints, there is a problem. I am also attaching the code here <https://github.com/DeepaMahm/misc/blob/master/cse_03_19_20.m>for your kind reference. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#41 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AB6CWOIMTBAKUN36NP7CNFDRKIE6XANCNFSM4LQMNEHQ> .