CVX Illegal Operation and Strange Zero Expression

1

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.

2

have you looked carefully at the result of the first iteration? Apparently, some iteration fails to solve, so the CXV variable values are populated with NaN, which results in invalid if used in a subsequent iteration. You can read some of my recent posts about the perils of unsafeguarded Successive Convex Approximation.

I haven’t looked to see what you’re doing with the Symbolic Toolbox, but be aware that you can’t use symbolic expressions in CVX code however you want - CVX doesn’t know how to deal with them.

3

Hi Mark,

So far I’ve not implemented any successive optimization procedure into the problem. I was planning to later and may have accidentally coded something like that in, but that was not my intention. What I meant by the iterations is that if I place a pause after the x within the loop it has no problem for say 30 or so iterations, but then an error is produced and NaN’s outputted.

The symbolic section is just to quickly generate the jacobian but then a matlabFunction is generated which can take regular values as inputs the same as any other user-produced matlab function. This may still be an issue but it certainly doesn’t cause any obvious issues.

4

The description in the 1st paragraph is too vague for me to understand.

I recommend you start with something simple, solve a instance of CVX, and carefully examine the results. When you have a complicated “mess”, it may be difficult to understand what’s going on, and why something fails.

5

I’m sorry Mark, I’ve solved the issue now. Denoting x as a variable rather than an expression fixed the issue. I am now getting a feasible output. The next step is to start including a cost function to minimize and path constraints such as normal force, heat rate and dynamic pressure. Thank you for your help.

6

Sorry for missing that. I didn’t look very carefully at your program.

7

That’s not a problem! I missed it all day today, but it’s all sorted now, thank you for your time.