Constrained Optimization Problem With Obj Function Which Is Numerical Solution To A System Of ODE's

The Problem:

Consider the following system of ODE's on an arbitrary time interval from [0,T]:

1. 
2. 
3. 
4. 

With Initial Conditions: , , ,

and where , , , , , ,

My goal is to minimize the final value of p(t) = p(T) by varying u and v subject to the following constraints:

The lower and upper bounds on the output of u(t) are:

The lower and upper bounds on the output of v(t) are:

The constraints are such that:

The sum of all outputs from u(t) may not exceed 300,( i.e. y = )

The sum of all outputs from v(t) may not exceed 10,( i.e. z = )

The Error:

While I'm not gettin an error message in the execution of the code, I keep getting answer structures which tell me the optimal solutions is one where u(t) = 0 and v(t) = 0 for all time. I suspect that the optimization problem solver might be considering p(T) as a constant and is recgonizes that no outputs of u(t) and v(t) will have an effect on a constant. If this is the case, how do I go about writing the problem so that the optimization solver consideres p(T) as a function of time and who's value depends on the dynamics of the ODE system above.

I'm also new to using matlab and have no prior coding experience, please try and explain things with that in mind.

My code is below:

The System of ODE's:

function dXdt = odeFunction(t,X)
 
 p = X(1);     
 q = X(2);      
 u = X(3);      
 v = X(4);      

 Theta = 0.1;          
 Gamma = 0.15;          
 b = 5.85;         
 d = 0.00873;       
 Mu = 0.02;        
 Zeta = 0.084;     
 Eta = 0;           

dpdt = -Zeta*p*log(p/q)-Theta*p*v;                          
dqdt = (b*p)-(Mu*q)-(d*(p)^(2/3)*q)-(Gamma*u*q)-(Eta*q*v);
dydt = u;
dzdt = v;

dXdt = [dpdt;dqdt;dydt;dzdt];
end
    

Solving The System With ode45

tRange = [0 T];
T = 10

p0 = 1;
q0 = 1;
u0 = 0;
v0 = 0;
IC = [p0;q0;u0;v0]; 

[tSol, XSol] = ode45(@odeFunction,tRange,IC);

p = real(XSol(:,1)); 
q = real(XSol(:,2)); 
y = real(XSol(:,3));
z = real(XSol(:,4));

pT = XSol(end,1); %This line gives the final output of p(t) at p(10)

Defining and Solving The Optimization Problem:

prob = optimproblem('Description','odeOptimization');

prob.Objective = pT

u = optimvar('u', 53, 'LowerBound', 0, 'UpperBound', 75);
v = optimvar('v', 53, 'LowerBound', 0, 'UpperBound', 2);

prob.Constraints.ymax = cumsum(u) <= 300;
prob.Constraints.zmax = cumsum(v) <= 10;

solve(prob)

[Matt J]

I suspect that the optimization problem solver might be considering p(T) as a constant and is recgonizes that no outputs of u(t) and v(t) will have an effect on a constant.

Yes, that is what is happening. If you run the code that computes pT, you will see that it is just a fixed number, 8.7923, with no dependence on u and v. And how could it depend on u and v? You do not even create u and v until well after pT is created.

T = 10;

tRange = [0 T];

p0 = 1;
q0 = 1;
u0 = 0;
v0 = 0;
IC = [p0;q0;u0;v0]; 

[tSol, XSol] = ode45(@odeFunction,tRange,IC);

p = real(XSol(:,1)); 
q = real(XSol(:,2)); 
y = real(XSol(:,3));
z = real(XSol(:,4));

pT = XSol(end,1) %This line gives the final output of p(t) at p(10)

pT = 8.7923

[Matt J]

It seems likely to me as well that you need constraints on p and q as well, since your ODEs involve the KL distance between p and q. Are they supposed to be non-negative distribution functions?

Answer 1

[Torsten]

You never use the 53 solution variables of u and v in your code.

And you solve the following ODEs for u and v

du/dt = u, u(0) = 0
dv/dt = v, v(0) = 0

which of course give u = v = 0 for all times as solution.

u and v don't have to be treated as differential variables, but as control variables.

Thus you have to solve the ODE problem

dpdt = -Zeta*p*log(p/q)-Theta*p*v;                          
dqdt = (b*p)-(Mu*q)-(d*(p)^(2/3)*q)-(Gamma*u*q)-(Eta*q*v);

with u and v given as

u = interp1(Tarray,Uarray,t)
v = interp1(Tarray,Varray,t)

where Tarray = linspace(0,T,53) and Uarray and Varray are your 53 solution variables u and v, respectively.

Answer 2

[Matt J]

Perhaps you might use fmincon as follows:

A=kron( eye(4) , ones(1,52)); b=[300;10;1;1]; %summation constraints
Aeq=[]; beq=[];

lb=zeros(52,4); ub=inf(size(lb)); 
ub(:,1)=75; % u-bounds
ub(:,2)=2; %v-bounds

lb(1,:)=[0,0,p0,q0]; %initial conditions
ub(1,:)=[0,0,p0,q0]; %initial conditions

objective=@(uvpq) uvpq(end,3);

uvpq0=rand(52,4); %You probably want a smarter initial guess than this.

uvpq=fmincon(objective, uvpq0, A,b,Aeq,beq,lb,ub,@nonlcon);

function [c,ceq]=nonlcon(uvpq)

     Theta = 0.1;          
     Gamma = 0.15;          
     b = 5.85;         
     d = 0.00873;       
     Mu = 0.02;        
     Zeta = 0.084;     
     Eta = 0;  
     
  c=[];

  vars=num2cell(uvpq,1);
  [u,v,p,q]=deal(vars{:});
  
  ceq(1)=-Zeta.*p.*log(p./q)-Theta.*p.*v - gradient(p);
  ceq(2)=(b.*p)-(Mu.*q)-(d.*(p).^(2/3)*q)-(Gamma.*u.*q)-(Eta.*q.*v) - gradient(q);
  
end