Solve optimization problem that iterates on parameters with constraints defined with numerical solutions to differential equations depending on those parameters

I want to search for an optimum in a problem where the constraints are numerical solutions of differential equations (these equations don't have analytical solutions) but the parameters I want to iterate on are parameters that appear in the differential equation so that the solver would need to:

• compute the numerical solution to the differential equations with current parameters
• check if my constraints on the numerical solutions are verified and compute
• if they are not, change parameters and compute the numerical solutions again to see if it gets closer
• once they are, try to minimize my objective function (which is also a function of these parameters) which to complicate even more is the maximum value in a vector

I don't have a lot of code to show because I don't know where to begin but basically I would like to do this :

objective: minimize (max(k))

constraints:

I=1;
k01=0.1;
theta_0=deg2rad(5);
k12=0.1;
theta_2=deg2rad(5);
Fm=10;
Fg=0.1;
l=30*10^(-3);
tspan = [0 5];
theta0 = [deg2rad(10) 0];

function dthetadt = odefcn(t,theta,k,I,theta_0,theta_2,Fm,Fg,l,n)
dthetadt = zeros(2,1);
dthetadt(1) = theta(2);
dthetadt(2) = (1/I)*(-k(1)*(theta(1)-theta_0)-k(2)*(theta(1)-theta_2)+Fm*l*sin(theta(1))-Fg*l*cos(theta(1)));
end

[t,theta] = ode45(@(t,theta) odefcn(t,theta,k), tspan, theta0);
theta(t_final)-theta_objective < tol

I guess come optimization toolboxes offer this possibility but I can't find one.

Thanks in advance.

Answer 1

[Alan Weiss]

Perhaps the example Fit an Ordinary Differential Equation (ODE) can help.

Alan Weiss

MATLAB mathematical toolbox documentation

[Johan Poccard]

Thank you for your answer.

I don't think the example applies to my case because I don't have data to fit the differenital equation solution to. I just want my solution y(t) equal to a certain value at a certain t. In other words I want:

[y1(tf),...,yn(tf)] = [yf1,...,yf] by iterating on para1,...,paraN

where the vector of solutions [y1(t),...,yn(t)] is defined by dd(yi)/dd(t)=f(y1,...,yn,ypara1,...,paraN) for i=1,...,n

[Alan Weiss]

Sorry, I don't understand your comment. If you fit the solution to those values, well, think of those values as the data and as far as I can tell the problems are identical. Obviously, I am missing something.

Alan Weiss

MATLAB mathematical toolbox documentation

[Johan Poccard]

Maybe I'm just not understanding the fit tool.

So using this fit tool, I could find the parameters that make my solutions go through one specific point(t,y) ?

I don't need a complete set of (t,y) to fit to ?

Can I set an objective funtion on the parameters for the fit if I have more parameters than points to go trough ? (I want the value of my parameters to be the lowest possible while still making the differential equation solution go trough specific points)

And is that a problem that I have a differential equation system as opposed to just one equation ?

[Alan Weiss]

If you fit your ODE to your data, there is no guarantee that the fit will be perfect, meaning that the fitted function exactly matches the data. But in any case, the fit will be the best possible, assuming that you reach the global minimum. You can, of course, use any model you like, with as many or few parameters as you want, though obviously having more parameters usually gives more flexibility for a better fit.

Please look at the example I linked. That is a system of three differential equations.

Alan Weiss

MATLAB mathematical toolbox documentation

[Johan Poccard]

Ok so I think I followed the example and adapted it accurately to my application but now I get the error 'LSQCURVEFIT requires all values returned by functions to be of data type double' which I don't understand since objective_fun returns pos which is similar to the one in the example.

Here's my code. I use an anonymous function objective_fun to pass extra parameters into the real objective function which is paramfun. paramfun itself calls dthetadt to comptute the diiferential equation system.

n=5;
k=[0.1 0.1 0.1 0.1 0.1 0.1];
tspan = [0 5];
theta0 = [deg2rad(linspace(5,-5,n)) zeros(1,n)];
Fm=10;

objective_fun=@(x,tspan)paramfun(x,tspan,n,theta0);

param0=[k, Fm];
xdata=5;
ydata=[0.4239    0.1086    0.0213   -0.1879   -0.3618];
[pbest] = lsqcurvefit(@(x,tspan)objective_fun,param0,xdata,ydata);

function dthetadt = odefcn(t,theta,k,Fm,I,Fg,l,n)
dthetadt = zeros(2*n,1);
clear vars i

dthetadt(1) = theta(n+1);
dthetadt(n+1) = (1/I)*(-k(1)*(theta(1))-k(2)*(theta(1)-theta(2))+Fm*l*sin(theta(1))-Fg*l*cos(theta(1)));

for i=2:n-1
    dthetadt(i) = theta(n+i);
    dthetadt(n+i) = (1/I)*(-k(i)*(theta(i)-theta(i-1))-k(i+1)*(theta(i)-theta(i+1))+Fm*l*sin(theta(i))-Fg*l*cos(theta(i)));
end

dthetadt(n) = theta(2*n);
dthetadt(2*n) = (1/I)*(-k(n)*(theta(n)-theta(n-1))-k(n+1)*(theta(n))+Fm*l*sin(theta(n))-Fg*l*cos(theta(n)));

end

function pos= paramfun(x,tspan,n,theta0)
k=x(1:n+1);
Fm=x(n+2);

Fg=0.1;
l=30*10^(-3);
I=1;

[~,pos]=ode45(@(t,theta) odefcn(t,theta,k,Fm,I,Fg,l,n),tspan,theta0);
end

[Alan Weiss]

When I run

tt = objective_fun(param0,tspan)

I get a value tt of size 53-by-10. The size of ydata is 1-by-5. Clearly, these are not the same. Take a look at the documentation of lsqcurvefit to see what your objective function should return, compared to ydata.

My guess is that you forgot that you have two values of time, 0 and 5, and you are giving data only for time 5.

Alan Weiss

MATLAB mathematical toolbox documentation

[Johan Poccard]

Thank you !