Nonlinear constraint on optimization using ODE - system state involved

Hello,

I need to incorporate a nonlinear constraint on my optimization problem but there is a ODE involved. The constraint involves a system state so that is why I am having some difficulty in incorporating it. I do not know if it can be done.

It is easy to write some constraint that depends only on the paramaters themselves or values that are available at t = 0 (e.g., initial conditions) but I do not know a priori the values of the states. So, adapting one MATLAB example from https://nl.mathworks.com/help/matlab/ref/ode23s.html let us say I am trying to optimize the two parameters in the vector p inside that ODE but some some reason I want to limit y(2) to a maximum value during optimization. Can someone help on this generic example?

tspan = [0 5];
y0 = [0 0.01];

function SSE = fobj(p,tspan,y0)

[t,y] = ode23s(@(t,y) odefcn(t,y,p), tspan, y0);
...  
    
SSE = sum((y - yexp).^2); % let us say we have some experimental data to define a sum of squared errors
    
end

function dydt = ode_example(t,y,p)

dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = p(1)*t.*y(1)+p(2);

end

function [c] = nlincon(y,p,max_value)

c = y(2)-max_value; % ?? how do I access the value of y(2) at all times outside the ODE function itself?

end

Answer 1

[Torsten]

I guess you have a time vector for which you have your data vector of measurements "yexp" ?

In nlincon, with these measurement times as tspan and with the actual parameter vector p, call the ODE function and define the vector c of size tspan as

c = y(:,2) - max_value

So defining tspan as [0 5] is not suitable. tspan should exactly contain the times at which "yexp" was measured.

[Gustavo Lunardon]

I am trying to contrain the values of the state y(2) at intermediates solution times, in a single shooting approach. Does your answer mean I have to solve the ode in the nonlinear constraint function call AND in the objective function call?

That would mean solving the same problem twice (computational effort might become an issue). Is there a workaround?

• If I make one file with both fobj and nonlincon functions would they be able to 'see' each other and avoid computing the same ode twice?

[Torsten]

I don't know the order and frequency in which the objective function and the constraint function are called.

But I doubt there is a rule you could exploit in order to save calls to the integrator.

At least, you can define an ODE function in which you solve the ODE. This function can then be called from the objective function and from the constraint function.

I don't understand how you want to constrain the values of the state y(2) in the shooting approach. Given the parameters p and the initial conditions, y(2) is given - it cannot be constraint. The only way is to force the optimizer to choose parameters that make y(2) bounded as required.

[Gustavo Lunardon]

That is exactly what I am trying to do, optimize the p vector to solve the problem under these additional constraints on the states. It is unfortunate I cannot avoid calling the integrator all the time in separate functions. That would save me a lot of simulation time.

It would be nice if I could make a function that gives

[obj,Aeq,beq,Aineq,bineq,lb,ub,nonlincon] = somefunction(inputs)

then fmincon or another solver would receive only that function instead of each piece of information in separate arguments

Maybe if I use the x = fmincon(problem) syntax then I could create a function that creates the problem to be optimized all at once?

[Torsten]

It is unfortunate I cannot avoid calling the integrator all the time in separate functions.

As said, you should call the integrator in one function that is called from the objective and the constraint function.

[Gustavo Lunardon]

Yes, I understood you suggested something of the form

function [t,y] = ode_solver(p,tspan,y0)
[t,y] = ode23s(@(t,y) odefcn(t,y,p), tspan, y0);
end

But I meant it was unfortunate I could not avoid solving the ode twice, even if I do like you said and call in integrator (ode_solver) from the objective and constraint function. Just to clarify. Thanks for the tip.

[Torsten]

Maybe MATLAB support can tell you more about order and frequency of calls to the input functions of fmincon in order to save computation time in your case.

I cannot.

[Torsten]

There was a similar question with the following answer given:

https://de.mathworks.com/help/optim/ug/objective-and-nonlinear-constraints-in-the-same-function.html

Hope it helps.