Question regarding passing input arguments in fmincon

Hi All,

I am solving an optimization problem.There are no equality constraints in my model. I want to pass ODEs as constraint to fmincon using,non-linear constraint, nlcon argument .

The function model contains odes

function dz = model(t, z, c)  
              dz(1) = ..
              dz(2) = ..

I'd like to ask how the ode intergrator has to be called in the place of nlcon in fmincon(@objective,p0,A,b,Aeq,beq,lb,ub,nlcon).

%p0 = Initial Value of Parameters
A = [];
b = [];
Aeq = [];
beq = [];
nlcon  = ode45(@(t,z)model(t,z,x), tSpan, z0); % Is this correct?
p = fmincon(@objective,p0,A,b,Aeq,beq,lb,ub,nlcon); 

Since the ode function is already called in nlcon , is it required to call the function model again inside the objective function?

function cost = objective(c,time_exp,expZ,tSpan,z0)
sol = ode45(@(t,z)model(t,z,c), tSpan, z0); % Is this step required?
ModelZ = deval(sol, time_exp);
cost = ModelZ-expZ;

Any help would be highly appreciated.

[Torsten]

Star Strider's answer should show you how to solve parameter fitting problems involving ODEs:

https://de.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting

[Deepa Maheshvare]

Thanks a lot for the response. I went through the details posted in the above-mentioned link.

Since the time series data was avilable in the example shown in the link, lsqcurvefit has been used.

In my case, I don't have the time series experimental data.

My real system is quite large, 100 odes , I am trying to perform global optimization.The data that I have from the experiments is steady state values of the variables z1,z2 and so on.

Probably my objective function should be modified as,

function cost = objective(c,expZ,tSpan,z0)
sol = ode45(@(t,z)model(t,z,c), tSpan, z0); % Is this step required?
ModelZ = sol(end,:)%  the solutions obtained from model at end point of tSpan
cost = ModelZ-expZ; %expZ contains steady state values of z1,z2.. 

since time_exp is the steady state time(not a vector).

Please note: I have upper and lower bounds (lb,ub)defined on the parameters that have to be estimated.

Could you please suggest how to proceed ?

[Torsten]

Use lsqcurvefit. Here you can also prescribe lower and upper bounds for the variables.It doesn't matter that ModelZ are steady-state values. Just proceed as in the given link.

[Deepa Maheshvare]

I'm actually trying to follow the description and example given in the post here.

I'd like to know whether the approach given in the above link wouldn't work.

[Deepa Maheshvare]

I had a look at this post. It explains a bit about the differences in fmincon and lsqcurvefit.

It states that "The simplest answer is that fmincon uses fancier algorithms than lsqnonlin and lsqcurvefit because fmincon must be able to deal with nonlinear constraints, whereas lsqnonlin/lsqcurvefit do not. "

In my case, I want to use the differential equations as nonlinear constraints.

My question is similar to the question posted here, where it has been advised to use noncoln option in fmincon to pass the odes as constraints in fmincon.

I'm really confused.

Any suggestions?

[Torsten]

If you want to restrict the values of the solution at all times to be in a certain range, not just the final values, you will have to use "fmincon" for your problem.

In "objective", return cost = sum((ModelZ-expZ).^2).

In "nlcon", choose an n-dimensional vector of time instants "tspan", solve the ODEs at these time instants and return a 2*n - dimensional constraint vector c with lb(i) <= c(2*i-1), ub(i) >= c(2*i) (i=1,...,n).

[Deepa Maheshvare]

Thank you very much. Yes, I want to do this "restrict the values of the solution at all times to be in a certain range, not just the final values".

May I ask a doubt regarding the implementation of this,"In "nlcon", choose an n-dimensional vector of time instants "tspan", solve the ODEs at these time instants "

tSpan = [0 2000] % in s 
nonlincon  = ode45(@(t,z)model(t,z,x), tSpan, z0);
% ode45(@(t,z)model(t,z,x), tSpan, z0); will return a matrix of dimension that is equal to(no. of time steps*m) m= the length of the vector z i.e. the number of variables z1,z2..,zm
x = fmincon(objective,x0,A,b,Aeq,beq,lb,ub,nonlincon);

From what you have explained, ncoln should be a vector (or matrix?)of dimension m*n i.e. 2*n (m = 2 , when z = z1,z2), where m is the number of variables and n is the number of timesteps. (Or 2 referes to upper bound and lower bound?).Please correct me if my understanding is wrong.

Could you please explain a bit more on this? 2*n - dimensional constraint vector c with lb(i) <= c(2*i-1), ub(i) >= c(2*i) (i=1,...,n).

The vector c, in my example in my original post, contains the initial values of the parametres that I want to estimate. I think here you are referring 'c' to the dimension of the vector returned

in nonlincon.

Also, in my original post, lb and ub is specified for the parameters in fmincon, for the paramters.

"lb(i) <= c(2*i-1), ub(i) >= c(2*i) (i=1,...,n)." Are these(lb, ub) for the variables(i.e z1,z2)?

[Torsten]

function [c,ceq] = nonlincon(x)
  tSpan = ...;
  z0 = ...;
  lb = ...;
  ub = ...;
  [T,Z] = ode45(@(t,z)model(t,z,x), tSpan, z0);
  n = numel(z0);
  m = numel(tSpan);
  for i = 1:n
    for j = 1:m
      c(2*((i-1)*m + j)-1) = Z(j,i) - ub(i);
      c(2*((i-1)*m + j))   = lb(i) - Z(j,i);
    end 
  end
end

If you set tSpan(end) to the value you also use in the objective function, you can call "fmincon" without setting lb and ub because "nonlincon" already accounts for these bounds at tSpan(end).

[Deepa Maheshvare]

Thanks a lot for helping me understand better.

I understand 'c' stores the difference between the values of the variable from the model and the maximum/minimum of the variable from the experimental observation.

c = Z(t,zi) - ub(i)

c = lb(i) - Z(t,zi)

This is computed for all variable at all time steps in the for loops.

Here lb and ub are the bounds of the variables.

I would like to know where the bounds of parameters have to be specified.

[Torsten]

Either where I wrote

lb = ...;
ub = ...;

or you can define all variables needed in "nonlincon" before the call to "fmincon" and pass them as

p = fmincon(@objective,p0,[],[],[],[],[],[],@(p)nonlincon(p,tSpan,z0,lb,ub));

[Deepa Maheshvare]

Apologies for the delay in my response,

What's confusing for me is

lb = ...;
ub = ...;

stores the bounds of the parameters.However, in the following loop

for i = 1:n
    for j = 1:m
      c(2*((i-1)*m + j)-1) = Z(j,i) - ub(i);
      c(2*((i-1)*m + j))   = lb(i) - Z(j,i);
    end 
  end

the difference is taken between the value variable,Z(j,i)(coputed from the ode solver) and the parameter,ub(i).

Could you please clarify this?

[Torsten]

The bounds lb and ub refered to in "nonlincon" do not store bounds for the parameters, but for the solution vector of your ODE (you wrote you wanted to restrict the values of the solution at all times to be in a certain range). If you want to use bounds for the parameters as well, you can do this directly in the call to "fmincon".

Note that lb, ub are of the same size as z0 in "nonlincon", not of the same size as the parameter vector x.

Note also that you should allocate c as

c = zeros(2*n*m,1);

before entering the double loop.

And set

ceq = [];

at the end of the function.

[Deepa Maheshvare]

Hi Torsten,

Thanks a ton for all the clarifications.

Since my model is in xml file , I am using SimBiology toolbox .I have created a model object using

modelObj = sbmlimport('Model.xml')

cs = getconfigset(modelObj);

cs.SolverType = 'ode15s';

cs.StopTime = 24;

[t, x, names] = sbiosimulate(modelObj);

There is a function call available for parameter estimation in SimBiology toolbox.

sbiofit(model,estimatedParam,'fmincon')

However, I am not sure how the input arguments of the estimation method,fmincon, has to be passed.I have a confusion here,

I've already specified p0,lb,ub in estimatedParam array with properties Name, InitialValue and Bounds

estimatedParams     = estimatedInfo(paramsToEstimate,'InitialValue',ParameterInitialValues,'Bounds',ParameterBounds)

Moreover, for the other arguments that have to be passed in the function call to fmincon,

fmincon(@objective,p0,A,b,Aeq,beq,lb,ub,nlcon)

p0 , lb,ub are already defined in estimatedParam in sbiofit(model,estimatedParam,'fmincon')

I'm not sure how to set up A=[],b=[],Aeq=[],beq = [] and define nlcon since the model is in xml format.

Could you please provide some suggestions on how to proceed?

[Torsten]

Sorry, but I don't have experience with the SimBiology toolbox.

Answer 1

[Alan Weiss]

I am not sure that I understand what you mean by "I want to pass ODEs as a constraint to fmincon." What about the ODE solution is a constraint? Do you want the minimum value of the solution to be above zero? Do you want the ODE solution to lie in a certain box? Once you explicitly define what you mean by the ODE solution is a "constraint," then I think that the answer to your question will be obvious. You will solve the ODE, get sol, and then write the constraint in terms of sol, maybe fun(deval(sol,tspan)) for some appropriate function fun and set of times tspan.

Alan Weiss

MATLAB mathematical toolbox documentation

[Deepa Maheshvare]

Explaining a bit about what I'm trying to do.

Simulating the model (without any optimization).

I have a system of stiff differential equations . I simulate the model using ode15s.The steady-state values that are generated from the integration solver isn't close to the experimental values. Therefore, I'm trying to estimate some parameters of my model while trying to

minimize the difference between the experimental values and values obtained from the model(solution from ode15s solver).

Do you want the ODE solution to lie in a certain box?

Yes, I have minimun and maxium values of z1 and z2 (observed from experiments).I want the solution that I obtain form the model to fall in this range(between min and max ).

[Alan Weiss]

I still do not understand precisely what you are trying to do, so cannot give you a very detailed answer.

If you want your steady-state values to be in a certain range, then you probably want to restrict deval(sol,lasttime), where lasttime is the final time. But if you just want to restrict the steady-state, then probably you shouldn't solve an ODE anyway, just the steady-state equations, which would be an fsolve call, finding the location where dz = zeros(size(z0)).

Maybe you want to restrict the values of the solution at all times to be in a certain range, not just the final values. In that case, you would set constraints on deval(sol,tspan), where tspan is the set of times that you want to restrict.

Maybe you are trying to fit your ODE solution to some observed (measured) values as a function of time. In that case, follow the procedures in this example.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Deepa Maheshvare]

Hi,

I'm sorry , I haven't explained it clearly.

"Maybe you are trying to fit your ODE solution to some observed (measured) values as a function of time."

No, the experimental measurements are all steady state values. For example , for the variable z1, there are 5 experimental studies that report the steady state value from experiment(z1_exp1,z1_exp2,z1_exp3,z1_exp4,z1_exp5). This gives me a range i.e. bound

z1_max and z1_min.

With z1_max and z1_min as the constraint , I wish to do the following

" restrict the values of the solution at all times to be in a certain range, not just the final values".Therefore, I am not sure whether it is appropriate to use lsqcurvefit ,used in the suggestions provided,for my problem.

[Alan Weiss]

Sorry, I am unable to understand what you are trying to do, so this will be my last reply. I suggest that you write out equations or inequalities describing your constraint. Seriously, write them down. Then, given those equations or inequalities, you can write a corresponding nonlinear constraint function in fmincon syntax, restricting the result by using deval on the solution at relevant times. If you don't understand what I am suggesting, sorry, maybe someone else can help.

Alan Weiss

MATLAB mathematical toolbox documentation