dsolve unable to convert to type 'branch'

I'm trying to find the optimal control of a system of differential equations subject to some cost function. The system of differential equations is given below:

The cost function is

The control is bang-bang (see below), and the function g is described below. The final time T is left free but we assume for the sake of example that . Model parameters include and γ. (which appears above) and (which appears below) are simply parameters denoting cost. The differential equations pertaining to the adjoint vector are

Where is the adjoint vector. The optimal control u is bang-bang and such that for and otherwise. Furthermore, the function g adjusts the cost associated with and . That is, for then , and for . is simply a scalar. Since I am trying to find the optimal control, this requires that the terminal (or end) conditions of the adjoint variables be . The initial conditions of the state variables are and . The end conditions of the state variables are left free. I tried solving this system with the following script (adapted from another optimal control question on this site):

clear all;
% helpful variables
beta = 2;
alpha = 0.1; 
gamma = 0.1;
M = 20;
IHat = 40;
T = 5;
k1 = 10;
k2 = 20;
c = 2;
N = 100;

% State equations
syms x1 x2 x3 p1 p2 p3 u;
Dx1 = -beta*x1*(x2+x3)/N + M*(x2/(x1+x2))*(p1 > (p2 + k1*(x1+x2)/x2)) + alpha*x3;
Dx2 = beta*x1*(x2+x3)/N - gamma*x2 - M*(x2/(x1+x2))*(p1 > (p2 + k1*(x1+x2)/x2));
Dx3 = gamma*x2 - alpha*x3;
% Cost function inside the integral
syms g;
g = -(k1*u + ((x2+x3) > IHat)*c*k2*(x2+x3) + ((x2 + x3) <= IHat)*k2*(x2+x3));
% Hamiltonian
syms p1 p2 p3 H;
H = g + p1*Dx1 + p2*Dx2 + p3*Dx3;
% Costate equations
Dp1 = -diff(H,x1);
Dp2 = -diff(H,x2);
Dp3 = -diff(H,x3);

% convert symbolic objects to strings for using 'dsolve'
eq1 = strcat('Dx1=',char(Dx1));
eq2 = strcat('Dx2=',char(Dx2));
eq3 = strcat('Dx3=',char(Dx3));
eq4 = strcat('Dp1=',char(Dp1));
eq5 = strcat('Dp2=',char(Dp2));
eq6 = strcat('Dp3=',char(Dp3));
sol_h = dsolve(eq1,eq2,eq3,eq4,eq5,eq6);
% use boundary conditions to determine the coefficients
% x1(0)=N-sum(I0); x2(0)=I0(1); x3(0)=I0(2);
% p1(T) = 0; p2(T) = 0; p3(T) = 0
conA1 = 'x1(0) = 90';
conA2 = 'x2(0) = 5';
conA3 = 'x3(0) = 5';
conA4 = 'p1(5) = 0';
conA5 = 'p2(5) = 0';
conA6 = 'p3(5) = 0';
sol_a = dsolve(eq1,eq2,eq3,eq4,eq5,eq6,conA1,conA2,conA3,conA4,conA5,conA6);

But I received the following error (when evaluating sol_h):

> In dsolve (line 121) 
Error using symengine
Unable to convert to type 'branch'.

Error in mupadengine/evalin_internal (line 113)
                res = mupadmex(statement,output_type{:});

Error in dsolve>mupadDsolve (line 327)
sys = [sys_sym reshape(evalin_internal(symengine, sys_str), 1, [])];

Error in dsolve (line 183)
sol = mupadDsolve(args, options);

I was unable to find another question with the same issue of conversion to type 'branch' so I figured I'd make a new question addressing this. Any help would be greatly appreciated.

[Walter Roberson]

Is there anything we can assume about the values? For example can we assume that the functions are completely real-valued ?

Answer 1

[Walter Roberson]

You had numerous mistakes in your program.

• You were trying to multiply the results of logical conditions by numeric values. In the symbolic toolbox, the result of logical tests is symtrue or symfalse, which are not convertable to numeric values. To get around this, I introduced a helper function symlogical which uses piecewise() to return 0 or 1
• You did not define any of your symbolic functions as being functions, but you differentiated them. Differentiating a symbol does not give the same result as differentiating a function, in ways that are important for differential equations
• you use variable names that are not obviously connected to your equations. Nothing I could do about that other than fix your syntax errors and trusting that you got the right mathematics (while being morally convinced that you did not)
• You differentiate H with respect to functions! That is not valid in normal Calculus, only in Calculus of Variations. You can only differentiate one function with respect to another when you can be sure that there are no non-obvious dependencies, which you can seldom be sure of when you are working with unknown functions. The closest you can get is to proceed as-if they are independent, and produce a hypothetical result, and then go back later and substitute and verify that indeed they are independent. I deliberately emit warning messages when the operation is done, as a reminder that the operation is suspicious
• you used character vectors for dsolve(). That is still tolerated in current releases, but expect it to go away soon. And with the other corrections needed... using the character vectors is likely to give the wrong results. Code has been adjusted to Don't Do That.
• When you run the resulting code, you will notice it says that it cannot find solutions. Sometimes that is just the way things are, there are a lot of differential equations that have no known explicit solution
• But in your case, the piecewise() that are being used is messing things up. dsolve() does not deal with piecewise() very well. What you need to do is go back and rewrite your comparisons that wrapped in symlogical() to instead be heaviside() calls: symlogical(A>B) -> heaviside(A - B) but pay attention to the value of heaviside(0) corresponding to A==B exactly, and see sympref('HeavisideAtOrigin')

% helpful variables
beta = 2;
alpha = 0.1; 
gamma = 0.1;
M = 20;
IHat = 40;
T = 5;
k1 = 10;
k2 = 20;
c = 2;
N = 100;

% State equations
syms x1(t) x2(t) x3(t) p1(t) p2(t) p3(t) H u g;
symlogical = @(CONDITION) piecewise(CONDITION, 1, 0);
Dx1 = -beta*x1*(x2+x3)/N + M*(x2/(x1+x2))*symlogical(p1 > (p2 + k1*(x1+x2)/x2)) + alpha*x3;
Dx2 = beta*x1*(x2+x3)/N - gamma*x2 - M*(x2/(x1+x2))*symlogical(p1 > (p2 + k1*(x1+x2)/x2));
Dx3 = gamma*x2 - alpha*x3;

% Cost function inside the integral
g = -(k1*u + symlogical((x2+x3) > IHat)*c*k2*(x2+x3) + symlogical((x2 + x3) <= IHat)*k2*(x2+x3));

% Hamiltonian
H = g + p1*Dx1 + p2*Dx2 + p3*Dx3;

% Costate equations
syms X1 X2 X3
Dp1 = -subs(diff(subs(H,x1,X1),X1),X1,x1);
warning('differentiating function H with respect to function x1 is not typically justifiable!')
Dp2 = -subs(diff(subs(H,x2,X2),X2),X2,x2);
warning('differentiating function H with respect to function x2 is not typically justifiable!')
Dp3 = -subs(diff(subs(H,x3,X3),X3),X3,x3);
warning('differentiating function H with respect to function x3 is not typically justifiable!')

% dsolve() does not like strings, and they do not give the right solutions.
eq1 = diff(x1) == Dx1;
eq2 = diff(x2) == Dx2;
eq3 = diff(x3) == Dx3;
eq4 = diff(p1) == Dp1;
eq5 = diff(p2) == Dp2;
eq6 = diff(p3) == Dp3;
eqns = [eq1,eq2,eq3,eq4,eq5,eq6];
sol_h = dsolve(eqns);

% use boundary conditions to determine the coefficients
% x1(0)=N-sum(I0); x2(0)=I0(1); x3(0)=I0(2);
% p1(T) = 0; p2(T) = 0; p3(T) = 0
conA1 = x1(0) == 90;
conA2 = x2(0) == 5;
conA3 = x3(0) == 5;
conA4 = p1(5) == 0;
conA5 = p2(5) == 0;
conA6 = p3(5) == 0;
conds = [conA1,conA2,conA3,conA4,conA5,conA6];
sol_a = dsolve([eqns, conds]);

[Cole Butler]

Thank you very much for your help, Walter. Even after making the change you suggested related to the heaviside function, it seems that no "explicit" solutions were found... Given that an explicit solution does not exist, is there a way of obtaining an implicit solution in MATLAB? I assume the format of such a solution would be an array of numbers corresponding to each variable at time t, but this is perfectly fine in my case.

[Walter Roberson]

If you follow the first example at odeFunction() you can produce anonymous functions to use with ode45()

However, your conditions appear to depend upon the evolution of the function, not just upon time. ode45() requires that the function be continuous for the time interval that it is to be integrated over. You will need to use event functions to detect when the conditions change, and tell it to exit integration at that point. And then you restart with the old boundary conditions and with any adjustments needed to the function so that it uses the new values. This is a bit of a nuisance...

[Cole Butler]

I see. My familiarity with MATLAB isn't that advanced to do it confidently, so it seems that an alternative method would be more appropriate. Thank you nonetheless for all your help.