How should i apply optimal control to my model? Following i have attached the code which i used to find maximum value of x3 in my model.

%%%%%%CODE TO FIND OPTIMUM (MAXIMUM) VALUE OF X3

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

function y = objfun(params)
options = odeset('InitialStep',0.0001,'RelTol',1e-09);
[~,Y] = ode15s(@(t,X)Scrpt1(t,X,params),[0 60],[9e-13,0],options);
y = -Y(end,2); % Negative because you want to maximize
% I assume that you want the value of Y at the last time
end

%%%%%%SOlUTION CODE %%%%%%

lb = [1e-9,4,0.1,4];
ub = [1e-6,15,2,10];
x0 = [1e-7,5,1,4];
opts = optimoptions("fmincon","PlotFcn","optimplotfval");
[sol,fval] = fmincon(@objfun,x0,[],[],[],[],lb,ub,[],opts)
% The true function value is -fval because you negated y

Answer 1

[Brahmadev]

Hi @Darshna Joshi,

As per my understanding, you would like to apply optimal control on the coefficients 'x1', 'x5', 'kst' and 'x0' that are provided by the 'params' parameter to the 'Scrpt1' function. For applying optimal control on these variables while maximizing the objective function, you can follow the following steps:

1. Define the control range of the parameters. See one example below:

% Assume we have 10 time intervals for control discretization
num_intervals = 10;

% Define the initial range for the control inputs at each interval. For
% example:
x1_control_sequence = linspace(1e-9, 1e-6, num_intervals);

2. Modify the objective function such that the 'ode15s' solver function is called iteratively with different time intervals. Change the example psudo-code below as per requirements.

% Defining time steps
dt = 60 / num_intervals;

% Initialize state and time
X = [9e-13, 0];
t = 0;

 % Simulate the system with piecewise constant control inputs
    for i = 1:num_intervals
        params = [x1_control_sequence(i); x5_control_sequence(i); kst_control_sequence(i); x0_control_sequence(i)];
        options = odeset('InitialStep',0.0001,'RelTol',1e-09);
        [T,Y] = ode15s(@(t,X)Scrpt1(t,X,params),[t t+dt],X,options);
        X = Y(end, :); % Update the state for the next interval
        t = t + dt; % Update time for the next interval
    end
    
    % Objective is the final value of x3 (maximization)
    y = -X(2);
        

3. Solve the Optimization Problem using the 'fmincon' function and setting the optimization options using 'optimoptions'.

Hope this helps in resolving your query.