How to set up parameter estimation in fmincon

Hi All,

I have the following dynamical systems. Equation 1 represents the exact dynamics of a system and equation 2 is the approximate dynamics that will give the same time course profiles as equation 1. To get the same time course profiles I have to determine D_hat in equation 2. I'm trying to solve this as a parameter estimation problem.

Dhat0 = %input vector 

% fun   = @objfun;

% [Dhat,fval] = fminunc(fun, Dhat0)


%% lsqnonlin

Dhat = lsqnonlin(@(Dhat) objfun(Dhat),Dhat0)



function f = objfun(Dhat)


%% Integrator settings

tspan = %tspan 

options = odeset('abstol', 1e-10, 'reltol', 1e-9);


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



%% generate approximate solution

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



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


%% objective function for lsqnonlin
    
f  = phi - phi_tilde

end

Using the above code, I could estimate D_hat in lsqnonlin.

Now, I am trying to solve this as an optimal conrol problem by using the equation (1) as dynamical constraints (non-linear equality constraints)

in fmincon.

I'm trying to set up the problem like the following

nonlcon = @defects;

Dhat= fmincon(@objfun,Dhat0,A,b,Aeq,beq,lb,ub,nonlcon)

For my system, A,b,Aeq,beq,lb,ub = []

But, I am not sure from where to pass these arguments for defects(dt,x,f).

function [c ceq] = defects(dt,x,f)
    
% ref: https://github.com/MatthewPeterKelly/OptimTraj
    
% This function computes the defects that are used to enforce the
    
% continuous dynamics of the system along the trajectory.
    
%
    
% 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
    
%
    
% OUTPUTS:
    
%   defects = [nState, nTime-1] = error in dynamics along the trajectory
    
%   defectsGrad = [nState, nTime-1, nDecVars] = gradient of defects



nTime = size(x,2);


idxLow = 1:(nTime-1);

idxUpp = 2:nTime;


xLow = x(:,idxLow);

xUpp = x(:,idxUpp);


fLow = f(:,idxLow);

fUpp = f(:,idxUpp);


% This is the key line:  (Trapazoid Rule)

defects = xUpp-xLow - 0.5*dt*(fLow+fUpp);

ceq = reshape(defects,numel(defects),1);

c = [];

end

end

Any suggestions will be really helpful

[Deepa Maheshvare]

Any suggestions?

Answer 1

[Alan Weiss]

I'm not sure, but I think what you are asking is how to pass several arguments as control variables, and maybe how to pass extra parameters (those that are not control variables, meaning you don't optimize over those variables).

If all of your variables dt, x, and f are control variables, then as documented you must put all of them in one array, typically called x, but I don't want to confuse it with your x variable, so I'll call the control variable v. Basically, you want v = [dt(:);x(:);f(:)]; In other words, v is a vector whose first component(s) are the dt entries, the next set are the x entries, and the last bits are the f entries. You code would look something like this:

function [c ceq] = defects(v)
    J = (length(v) - 1)/2;
% INPUTS:
%   dt = time step (scalar)
    dt = v(1);
%   x = [nState, nTime] = state at each grid-point along the trajectory
    x = v(2:(J+1));
%   f = [nState, nTime] = dynamics of the state along the trajectory
    f = v((J+1):end);
% OUTPUTS:
%   defects = [nState, nTime-1] = error in dynamics along the trajectory    
%   defectsGrad = [nState, nTime-1, nDecVars] = gradient of defects

nTime = size(x,2);

idxLow = 1:(nTime-1);
idxUpp = 2:nTime;

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

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

% This is the key line:  (Trapazoid Rule)
defects = xUpp-xLow - 0.5*dt*(fLow+fUpp);

ceq = reshape(defects,numel(defects),1);
c = [];

end

end

If I am wrong and you don't want to optimize all of these variables, then see Passing Extra Parameters.

Alan Weiss

MATLAB mathematical toolbox documentation

[Deepa Maheshvare]

Hi Alan,

I am not sure if this will work,

` In other words, v is a vector whose first component(s) are the dt entries, the next set are the x entries, and the last bits are the f entries.`

The defects function that I had posted in my original post is from here https://github.com/MatthewPeterKelly/OptimTraj/blob/master/trapezoid.m (please refer to function named computeDefects). I'd mentioned this reference as a comment in the code, but didn't explicitly state it. Sorry about that.

% 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

[nState, nTime] and [nState, nTime] refers to the size of the varible x and f. From what I understand, these are matrices.

In the link mentioned above, I am not sure how to input arguments are passed for the function computeDefects (line 113) defined in the file trapezoid.m .

I think x and f are computed by solving the system dynamics that are passed as nonlinear quality contraints, obtained by trapezoidal collocation) to fmincon(in my example, eq(1) defined in the original post)

In line 71 of trapezoid.m ,

 % Trapazoid integration calculation of defects:
    problem.func.defectCst = @computeDefects;

computeDefects is defined assigned to problem.func.defectCst. And problem is a struct passed to fminconfmincon(problem), https://in.mathworks.com/help/optim/ug/fmincon.html.

My confusion is how are the input arguments passed (from which function of optimtraj or is ti done internally in fmincon) to @computeDefects defined as function defects in my original post.

Please let me know if I am still not clear

[Alan Weiss]

You need to ensure that the arrays are passed properly once extracted from the v vector. You might need to reshape things internally:

function [c ceq] = defects(v)
    J = (length(v) - 1)/2;
% INPUTS:
%   dt = time step (scalar)
    dt = v(1);
%   x = [nState, nTime] = state at each grid-point along the trajectory
    x = v(2:(J+1));
%   f = [nState, nTime] = dynamics of the state along the trajectory
    f = v((J+1):end);
% OUTPUTS:
%   defects = [nState, nTime-1] = error in dynamics along the trajectory    
%   defectsGrad = [nState, nTime-1, nDecVars] = gradient of defects
nState = 5; % PUT IN THE CORRECT VALUE HERE
nTime = 12; % PUT IN THE CORRECT VALUE HERE
x = reshape(x,nState,nTime);
f = reshape(f,nState,nTime);
% Proceed with your calculations below

If I misunderstand you, I am sorry, but at some point you have to understand what sizes the arrays are and build that understanding into your function.

Alan Weiss

MATLAB mathematical toolbox documentation

[Deepa Maheshvare]

Hi Alan,

`If all of your variables dt, x, and f are control variables, then as documented you must put all of them in one array, typically called x, but I don't want to confuse it with your x variable, so I'll call the control variable v. `

Actually, dt, x, and f are not control variables. I've also attached an image for your kind reference. y dot = f[y,u,t]. In my case, dx/dt = f . (phi is replaced in place of x in equation1 presented in my original post).

If I am still not clear, I will add a detailed documentation.

[Alan Weiss]

You must be able to determine exactly what your control variables are. If those variables you mentioned are not control variables, then you should not be using them in your nonlinear constraint function except possibly as extra parameters (as in Passing Extra Parameters).

Please do not attach detailed documentation. Just say what are the control variables. If you cannot explain concisely exactly what are your control variables, then I will not be able to help.

Alan Weiss

MATLAB mathematical toolbox documentation

[Deepa Maheshvare]

Hi Alan,

My control variables are defined in equation2 in original post (re-written in terms of x below).

Matrix, M is known.I want to determine that minimizes the cost function.

The non-linear constraints have been written by discretizing equation (2) using trapezoidal collocation.

After what you suggested, I have defined function [c ceq] = defects(v) that accepts as input argument v and the x and dx values necessary for setting up the constraints are obtained by calling a nested function defined within defects` .This nested function has accepts input arguments ode15s(@(t,phi_tilde) fun(t,x, )`.

I hope this is right

[Deepa Maheshvare]

However, I ran into a dimension error and I have opened a new thread here