Optimization using fmincon ODE Matlab

Question

The problem is to find the optimum(maximum) value of x3 in range of (-8e-4 to 2e-4) by varying kst,x1,x5 and xo)

x5=5    %Input 2 (Input 2 is a state variable and could vary in range of 4 to 15 while performing 
                  optimization)
kst=1   %Input 3 (Input 3 is in terms of rate constant, it could vary from 0.1 to 2)
xo=4    %Input 4 (Input 4 is a state variable and could vary in range of 4 to 10)
x1=1e-7 %Input 1 could vary from 1e-9 to 1e-6

1. Script file

    function rest = Scrpt1(t,X)
    x2 = X(1); 
    x3 = X(2); 
    
    %Parameters
    
    if t<15
    
    x1 = 1e-7; %Input 1 could vary from 1e-9 to 1e-6
    
    else x1 = 0;
    
    end
    
    x5=5 %Input 2 (Input 2 is a state variable and could vary in range of 4 to 15 while performing optimization)
    
    kst=1 %Input 3 (Input 3 is in terms of rate constant, it could vary from 0.1 to 2)

    xo=4 %Input 4 (Input 4 is a state variable and could vary in range of 4 to 10)
    
    k1 = 6e7;

    km1 = 0.20;

    km4 = 0.003;

    k3 = 2500.00;

    k4 = km4/9;

    km3 = km1;

    LAP=1.5

% Differential equations

    dx2dt = km1*x3 + km3*LAP - k1*x1*x2 + km4*x3 - k4*x2;

    dx3dt = k1*x1*x2 - km1*(x3+x5+xo) - k3*x3*kst;

    rest = [dx2dt; dx3dt];

    end

2. Function file for ODE solution

    options = odeset('InitialStep',0.0001,'RelTol',1e-09);
    
    [T,Y]=ode15s(@Scrpt1,[0 60],[9e-13,0],options);
    
    X3= Y(:,2);
    
    plot(T,X3)

How to use fmincon or any other optimization solver for this to solve the mentioned optimization problem of finding maximum value of x3. For which values of x5,kst,xo,x1 we get maximum x3?

Answer 1

First you must add the values you want wo optimie as parameters of your coupled diffrential equations:

function rest = Scrpt1(t,X,X_opt)
    
    x5=X_opt(1);
    kst=X_opt(2);
    xo=X_opt(3);
    x1=X_opt(4);  
    x2 = X(1); 
    x3 = X(2); 
    
    %Parameters
    
    if t>=15
    x1 = 0;  
    end
    
    k1 = 6e7;
    km1 = 0.20;
    km4 = 0.003;
    k3 = 2500.00;
    k4 = km4/9;
    km3 = km1;
    LAP=1.5;

% Differential equations

    dx2dt = km1*x3 + km3*LAP - k1*x1*x2 + km4*x3 - k4*x2;
    dx3dt = k1*x1*x2 - km1*(x3+x5+xo) - k3*x3*kst;
    rest = [dx2dt; dx3dt];
    end

then you have to wirte a wrapper function you want to minimize. Because you want to maxmize x3 you have to add an minus to your objective value.

function max_X3=fun(X_opt) 
    tspan=[0 60];
    y0=[9e-13,0];
    options = odeset('InitialStep',0.0001,'RelTol',1e-09);
    
    [~,y] = ode15s(@(t,y) Scrpt1(t,y,X_opt), tspan, y0,options);
    
    max_X3=-max(y(:,2));
end

Finally you can use fmincon like this:

% x5, kst, xo, x1
initial_search_point=[5, 1, 4, 1e-7]
lower_bounds=[4, 0.1, 4, 1e-9]
upper_bounds=[15, 2, 10, 1e-6]

fmincon(@fun,initial_search_point,[],[],[],[], lower_bounds,upper_bounds)