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

differential equationsnumerical integrationoptimization

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.

Best Answer

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

Alan Weiss

MATLAB mathematical toolbox documentation

Related Solutions

MATLAB: ‘LSQCURVEFIT requires all values returned by functions to be of data type double’ error

It took a bit of effort to figure out what you are doing, and get your code to run correctly, although it still needs for you to make it work as you want it to. I cannot do that. Only you ccan make those decisions.

Specifically, in ‘paramfun’ you need to decide what you want it to return, since it must be a vector the same size as ‘ydata’. As wrritten, the ode45 call returns a (5x10) matrix, and there is no way that is going to fit a (1x5) vector.

The revised code:

n=5;
k=[0.1 0.1 0.1 0.1 0.1 0.1];
ydata=[0.4239    0.1086    0.0213   -0.1879   -0.3618];
tspan = linspace(0, 5, numel(ydata));                                       % Since ‘tspan’ Is The Independent Variable, It Needs to Be The Same Size
theta0 = [deg2rad(linspace(5,-5,n)) zeros(1,n)];
Fm=10;
objective_fun=@(x,tspan)paramfun(x,tspan,n,theta0);
param0=[k, Fm];
[pbest] = lsqcurvefit(objective_fun,param0,tspan,ydata)                     % The ‘objective_fun’ Function Already Exists To Work With ‘lsqnonlin’ Correctly, ‘tspan’ Is ‘xdata’
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 posout= 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);
posout = pos(:,1).';                                                        % Choose A Column, Then Transpose It So That ‘posout’ & ‘ydata’ Are The Same Size
end

This has the virtue of running without error (although with one warning and a note that it stopped because lsqcurvefit exceeded the funciton evaluation limit, not that it had found a global minimum and a reasonably correct solution.

.

MATLAB: Solving non linear equation using fsolve

You are not passing theta_n and other parameters to root2d(). One fix is to make root2d a Nested Function. Also, you must choose a different initial point. phi0 = [10,10] worked for me.

function solve_it()
    beta_a = 10;
    alpha_n = 10 ;
    eta= 18;
    i = 18;
    theta_i = 10;
    theta_n = 10 ;
    
    fun = @root2d;
    phi0 = [10,10];
    phi = fsolve(fun,phi0);
    phi_n = phi(1),
    phi_i = phi(2),
    function F = root2d(phi)
    F = zeros(2,1);
    F(1) = (pi*sin((pi*(phi(1) + theta_n))/180)*cos((pi*phi(2))/180)*(cos((pi*theta_n)/180) + tan((pi*i)/180)*tan((pi*theta_i)/180)))/(180*sin((pi*phi(1))/180)*(cos((pi*(phi(1) + theta_n))/180)*cos((phi(2)*pi)/180) + sin((phi(2)*pi)/180)*tan((theta_i*pi)/180))^2) - (pi*cos((pi*phi(1))/180)*(cos((pi*theta_n)/180) + tan((pi*i)/180)*tan((pi*theta_i)/180)))/(180*sin((phi(1)*pi)/180)^2*(cos((pi*(phi(1) + theta_n))/180)*cos((pi*phi(2))/180) + sin((pi*phi(2))/180)*tan((pi*theta_i)/180)));
    F(2) = ((cos((pi*theta_n)/180) + tan((pi*i)/180)*tan((pi*theta_i)/180))*((pi*cos((pi*(phi(1) + theta_n))/180)*sin((pi*phi(2))/180))/180 - (pi*cos((pi*phi(2))/180)*tan((pi*theta_i)/180))/180))/(sin((pi*phi(1))/180)*(cos((pi*(phi(1) + theta_n))/180)*cos((phi(2)*pi)/180) + sin((phi(2)*pi)/180)*tan((theta_i*pi)/180))^2);
    end
end

Related Question