MATLAB: Dsolve unable to convert to type ‘branch’

differential equationsMATLABoptimization

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.

Best Answer

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]);

Related Solutions

MATLAB: Concerning some error For Loop in the optimal control problem

Inside a quoted string 'i' the symbolic engine will always interpret 'i' as sqrt(-1) and never as the current value of i

You should assume that you are no longer permitted to pass character vectors to dsolve and rewrite everything to symbolic expressions

You are doing some invalid calculations. Dx1 only makes sense if x1 is a function, but if it is a function then you cannot take the derivative of H with respect to it. Taking the derivative of a function with respect to a function is an invalid operation in normal calculus and requires Calculus of Variations instead. Unknown functions cannot be assumed to be independent, especially in a situation like this where the functions are being defined by their relationships to each other.

MATLAB: How to solve in MATLAB 2018b ??

Hi,

the following runs for me in 2018b:

clear all
% State equations
syms x1 x2 p1 p2 u;
Dx1 = x2;
Dx2 = -x2 + u;
% Cost function inside the integral
syms g;
g = 0.5*u^2;
% Hamiltonian
syms p1 p2 H;
H = g + p1*Dx1 + p2*Dx2;
% Costate equations
Dp1 = -diff(H,x1);
Dp2 = -diff(H,x2);
% solve for control u
du = diff(H,u);
sol_u = solve(du,u);
% Substitute u to state equations
Dx2 = subs(Dx2,u,sol_u);
% convert symbolic objects to strings for using 'dsolve'
eq1 = strcat('Dx1=',char(Dx1));
eq2 = strcat('Dx2=',char(Dx2));
eq3 = strcat('Dp1=',char(Dp1));
eq4 = strcat('Dp2=',char(Dp2));
sol_h = dsolve(eq1,eq2,eq3,eq4);
% use boundary conditions to determine the coefficients
%    case a: (a) x1(0)=x2(0)=0; x1(2) = 5; x2(2) = 2;
conA1 = 'x1(0) = 0';
conA2 = 'x2(0) = 0';
conA3 = 'x1(2) = 5';
conA4 = 'x2(2) = 2';
sol_a = dsolve(eq1,eq2,eq3,eq4,conA1,conA2,conA3,conA4);
% plot both solutions
figure(1);
ezplot(sol_a.x1,[0 2]); hold on;
ezplot(sol_a.x2,[0 2]);
ezplot(-sol_a.p2,[0 2]);    % plot the control: u=-p2

axis([0 2 -1.6 7]);
text(0.6,0.5,'x_1(t)');
text(0.4,2.5,'x_2(t)');
text(1.6,0.5,'u(t)');
xlabel('time');
ylabel('states');
title('Solutions comparison (case a)');
%    case b: (a) x1(0)=x2(0)=0; p1(2) = x1(2) - 5; p2(2) = x2(2) -2;
eq1b = subs(sol_h.x1,'t',0);
eq2b = subs(sol_h.x2,'t',0); 
eq3b = subs(sol_h.p1,'t',2) == subs(sol_h.x1,'t',2)-5;
eq4b = subs(sol_h.p2,'t',2) == subs(sol_h.x2,'t',2)-2;
sol_b = solve(eq1b,eq2b,eq3b,eq4b);
C1 = double(sol_b.C1);
C2 = double(sol_b.C2);
C3 = double(sol_b.C3);
C4 = double(sol_b.C4);
sol_b2 = struct('x1',{subs(sol_h.x1)},'x2',{subs(sol_h.x2)}, ...
                'p1',{subs(sol_h.p1)},'p2',{subs(sol_h.p2)});
figure(2);
ezplot(sol_b2.x1,[0 2]); hold on;
ezplot(sol_b2.x2,[0 2]);
ezplot(-sol_b2.p2,[0 2]);    % plot the control: u=-p2
axis([0 2 -0.5 3]);
text(0.9,0.5,'x_1(t)');
text(0.4,1,'x_2(t)');
text(0.2,2.5,'u(t)');
xlabel('time');
ylabel('states');
title('Solutions comparison (case b)');

Best regards

Stephan

Best Answer

Related Solutions

MATLAB: Concerning some error For Loop in the optimal control problem

MATLAB: How to solve in MATLAB 2018b ??

Related Question