I’m getting the following error message: Illegal operation: {constant} + {Invalid} in my CVX program. It seems to iterate through initially but eventually become an issue. As part of the code I’m outputting the Jacobian, B*u, State update © and unbounded controller EPS. The program has marked the Jacobian as a “cvx zero expression” which doesn’t sound right, and the others as affine which does seem correct. Here is the code.

```
%% CVX PREP + REF TRAJ
R0 = 6378e3; % Radius of Earth in m.
g0 = 9.81; % Acceleration due to gravity in m/s^2.
kQ = 1.65e-4*sqrt(g0*R0)^(3.15); % Heat transfer constant.
rho0 = 1.225; % Density of air in kg/m^3;
hs = 7200; % Reference altitude in m.
Omega = 7.2921150e-5; % Rotational rate of Earth.
no = 1;
xi = [100e3;0;0;7450;-0.5*pi/180;0;0];
xf = [25e3;12*pi/180;70*pi/180;0;-10*pi/180;90*pi/180;0];
r = xi(1);
th = xi(2);
ph = xi(3);
V = xi(4);
ga = xi(5);
ps = xi(6);
si = xi(7);
B = [zeros(6,1);1];
la = 10^8;
ti = 0; % Initial time in s.
tf = 42; % Final time in s.
N = tf*20; % Time segments.
dt = (tf - ti)/(N-1); % t delta.
xdim = numel(xi);
udim = 1;
ttest = ti:dt:tf;
switch no
case 1
rn = r/R0;
Vn = V/sqrt(R0*g0);
tin = ti/sqrt(R0/g0);
tfn = tf/sqrt(R0/g0);
dtn = (tfn-tin)/(N-1);
Omegan = Omega/sqrt(g0/R0);
case 0
rn = r;
Vn = V;
tin = ti;
tfn = tf;
dtn = (tfn-tin)/(N-1);
Omegan = Omega;
end
xin = xi./[R0,1,1,sqrt(R0*g0),1,1,1]';
xfn = xf./[R0,1,1,sqrt(R0*g0),1,1,1]';
m = 104305; % Mass in KG.
Ar = 391.22; % Surface area in m^2.
rx(:,1) = [rn;th;ph;Vn;ga;ps;si];
rx(:,2:N) = zeros(7,N-1);
for i = 1:N
xd = StateFunction(rx(1,i),rx(2,i),rx(3,i),rx(4,i),rx(5,i),rx(6,i),rx(7,i),no);
rx(:,i+1) = rx(:,i) + dtn*xd;
% if rx(1,i+1) < 0
% rx(1,i+1) = 0;
% break
% end
% if rx(4,i+1) < 0
% rx(4,i+1) = 0;
% break
% end
end
%% Constraints
sigmadotmax = 10*pi/180;
Qdmax = 1500;
qmax = 18000;
nmax = 2.5*g0;
sigmamax = 90*pi/180;
%% CVX
k = 1;
cvx_begin
% Variables
variables u(udim,N) ep(xdim,N)
expressions x(xdim,N+1)
% Cost Function
% minimize(-x(4,end) + la*sum(sum(ep)));
% Initial and Final conditions.
x(:,1) == xin;
x(1,N) == xfn(1); % Length from centre of Earth.
x(2,N) == xfn(2); % Longitude
x(3,N) == xfn(3); % Latitude.
x(5,N) == xfn(5); % Flight path angle.
x(6,N) == xfn(6); % Heading angle.
% Path constraints.
for i = 1:N
x(1,i) >= 0;
x(4,i) >= 0;
x(7,i) <= sigmamax;
x(7,i) >= -sigmamax;
u(:,i) <= sigmadotmax;
u(:,i) >= -sigmadotmax;
end
% Equations of motion.
for i = 1:N
x(:,i+1) == x(:,i) + dtn*(Afunc(rx(1,i),rx(2,i),rx(3,i),rx(4,i),rx(5,i),rx(6,i),rx(7,i))*x(:,i) + B*u(:,i) + C(rx(:,i),no) + ep(:,i));
a = dtn*Afunc(rx(1,i),rx(2,i),rx(3,i),rx(4,i),rx(5,i),rx(6,i),rx(7,i))*x(:,i)
b = dtn*B*u(:,i)
c = dtn*C(rx(:,i),no)
d = dtn*ep(:,i)
pause
clc
end
cvx_end
function [output] = C(x,no)
output = StateFunction(x(1),x(2),x(3),x(4),x(5),x(6),x(7),no)-Afunc(x(1),x(2),x(3),x(4),x(5),x(6),x(7))*x;
end
function f = StateFunction(varargin)
switch nargin
case 8
r = varargin{1};
th = varargin{2};
ph = varargin{3};
V = varargin{4};
gam = varargin{5};
ps = varargin{6};
sig = varargin{7};
no = varargin{8};
R0 = 6378e3;
g0 = 9.8;
hs = 7200;
Aref = 391.22;
switch no
case 1
Omega = 7.2921150e-5/sqrt(g0/R0);
case 0
Omega = 7.2921150e-5;
end
switch no
case 1
h = r*R0;
case 0
h = r;
end
rho0 = 1.225;
rho = rho0*exp(-h/hs);
m = 104305;
M = (V-4570)^2/(340^2);
alph = 40*(V>4570) + (40 - (0.20705*M))*(V<=4570);
CL = -0.041065 + 0.016292*alph + 0.00026024*alph^2;
CD = 0.080505 - 0.03026*CL + 0.86495*CL^2;
otherwise
syms V th gam ps r ph sig alph real
R0 = 6378e3;
g0 = 9.8;
hs = 7200;
Aref = 391.22;
no = 0;
switch no
case 1
Omega = 7.2921150e-5/sqrt(g0/R0);
case 0
Omega = 7.2921150e-5;
end
switch no
case 1
h = r*R0;
case 0
h = r;
end
rho0 = 1.225;
rho = rho0*exp(-h/hs);
m = 104305;
CL = -0.041065 + 0.016292*alph + 0.00026024*alph^2;
CD = 0.080505 - 0.03026*CL + 0.86495*CL^2;
end
switch no
case 1
L = (R0*rho*V^2*Aref*CL)/(2*m);
D = (R0*rho*V^2*Aref*CD)/(2*m);
Ln = L/g0;
Dn = D/g0;
case 0
Ln = (rho*V^2*Aref*CL)/(2*m);
Dn = (rho*V^2*Aref*CD)/(2*m);
end
f = [V*sin(gam);...
(V*cos(gam)*sin(ps))/(r*cos(ph));...
(V*cos(gam)*cos(ps))/r;...
-Dn-(sin(gam))/(r^2);...
((Ln*cos(sig))/(V)) + ((V^2-(1/r))*cos(gam))/(V*r);...
((Ln*sin(sig))/(V*cos(gam))) + (V*cos(gam)*sin(ps)*tan(ph))/r;...
0];
end
function A = Afunc(r,th,ph,V,gam,ps,sig)
%AF
% A = AF(R,TH,PH,V,GAM,PS,SIG)
% This function was generated by the Symbolic Math Toolbox version 8.7.
% 30-Apr-2021 09:41:31
M = (V-4570)^2/(340^2);
alph = 40*(V>4570) + (40 - (0.20705*M))*(V<=4570);
t2 = cos(gam);
t3 = cos(ph);
t4 = cos(ps);
t5 = cos(sig);
t6 = sin(gam);
t7 = sin(ps);
t8 = tan(ph);
t9 = sin(sig);
t10 = V.^2;
t11 = alph.^2;
t12 = 1.0./V;
t13 = 1.0./r;
t19 = r./7.2e+3;
t23 = alph.*1.6292e-2;
t25 = alph.*4.9299592e-4;
t14 = t13.^2;
t15 = t13.^3;
t16 = 1.0./t2;
t17 = 1.0./t3;
t18 = -t13;
t20 = -t19;
t24 = t11.*2.6024e-4;
t26 = t11.*7.8748624e-6;
t21 = exp(t20);
t22 = t10+t18;
t27 = t23+t24-4.1065e-2;
t28 = t27.^2;
t29 = t28.*8.6495e-1;
t30 = -t29;
t31 = t25+t26+t30-8.174762689999999e-2;
mt1 = [0.0,-V.*t2.*t7.*t14.*t17,-V.*t2.*t4.*t14,t6.*t15.*2.0-t10.*t21.*t31.*3.190726049140075e-7,t2.*t12.*t15-V.*t5.*t21.*t27.*3.190726049140075e-7-t2.*t12.*t14.*t22,-V.*t2.*t7.*t8.*t14-V.*t9.*t16.*t21.*t27.*3.190726049140075e-7,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,V.*t2.*t7.*t13.*t17.^2.*sin(ph),0.0,0.0,0.0,V.*t2.*t7.*t13.*(t8.^2+1.0),0.0,t6,t2.*t7.*t13.*t17,t2.*t4.*t13,V.*t21.*t31.*4.594645510761708e-3];
mt2 = [t2.*t13.*2.0+t5.*t21.*t27.*2.297322755380854e-3+(t2.*t18.*t22)./t10,t2.*t7.*t8.*t13+t9.*t16.*t21.*t27.*2.297322755380854e-3,0.0,V.*t2,V.*t6.*t7.*t17.*t18,V.*t4.*t6.*t18,-t2.*t14,t6.*t12.*t18.*t22,V.*t6.*t7.*t8.*t18+V.*t6.*t9.*t16.^2.*t21.*t27.*2.297322755380854e-3,0.0,0.0,V.*t2.*t4.*t13.*t17,V.*t2.*t7.*t18,0.0,0.0,V.*t2.*t4.*t8.*t13,0.0,0.0,0.0,0.0,0.0,V.*t9.*t21.*t27.*(-2.297322755380854e-3),V.*t5.*t16.*t21.*t27.*2.297322755380854e-3,0.0];
A = reshape([mt1,mt2],7,7);
end
```

Does anyone know what’s going wrong? If you need more information please let me know.