Variables and objective value zero:how to analyze the results from CVX (solver sdpt3)


#1

Hello all,

I would really appreciate if someone can please point out the errors in my code. My variables and objective values are coming out to be zero and nearly zeros. The problem I am trying to reproduce is a spacecraft reorientation problem (https://ieeexplore.ieee.org/document/1271671/). The problem and my code are as follows:

clc;clear all;close all;

ti=0;tf=20;dt=1;t=[0:dt:tf]’;
n =((tf-ti)/dt);%divide t0-tf into n intervals

J=diag([100,200,300]);

x=zeros(n,10);
u=x(2:end,1:3);
w=x(:,4:6);
q=x(:,7:10);
X=x(2:end,:);

% Initial and final conditions of angular velocity and quaternions
q(1,:)=[0.5, 0.5, 0.5, 0.5];%q(1,:)=q(1,:)/norm(q(1,:));
w(1,:)=[0.2,0.2,0.2]*pi/180;%[0.2,0.2,0.2]pi/180;[sin(0.01) sin(0.0085) cos(0.0085)]
wdot(1,:)=zeros(3,1)’;% calculated analytically from w(1,: )
u(1,:)= (J
wdot(1,:)’+crossm(w(1,:)’)Jw(1,:)’)’;%zeros(1,3);% calculated from initial values of q and w

q(end,:)=[0.0258, 0.0258, 0.999, 0.0258];%q(end,:)=q(end,:)/norm(q(end,:));
w(end,:)=[0,0,0];

x(1,:)=[ zeros(1,3), w(1,:), q(1,:)];
x(end,:)=[u(end,: ), w(end,: ), q(end,:)];

G1=[eye(3,3) zeros(3,7)];%dim=3x10
gamma1=03;% control bound (rad/s^2)

G2=[zeros(3,3) eye(3) zeros(3,4)];%dim=3x10
gamma2=.3;% angular velocity bound (rad/s)

xs_I=[0; 0; 1];% inertial coordinates of sun vector
yc_B=[0.750;0.4330;0.500];%spacecraft camera vector in body coordinates

A =xs_Iyc_B’+yc_Bxs_I’-(xs_I’yc_B+cos(50pi/180))*eye(3);%dim=3x3
b = [0 -xs_I(3) xs_I(2);xs_I(3) 0 -xs_I(1);-xs_I(2) xs_I(1) 0]yc_B;%dim=3x1
d = xs_I’yc_B-cos(50pi/180);%dim=1x1
Atilda = [A b;
b’ d];
% testAtilda_sym = issymmetric(Atilda);
% testAtilda_psd = isPositiveDefinite(Atilda);
eig_atilda=eig(Atilda);
tol=1e-12;
mu=max(tol,abs(min(eig(Atilda))));%abs(min(eig(Atilda)))
Atilda1=Atilda+mu
eye(4);% making the matrix in non-convex constraint positive semi-definite
eig_atilda1=eig(Atilda1);

H = [zeros(6,6) zeros(6,4);zeros(4,6) Atilda1];%dim=10x10
C=[eye(3) zeros(3,1)];

Q_n=[q(end, 4) q(end,3) -q(end,2) -q(end,1);
-q(end,3) q(end,4) q(end,1) -q(end,2);
q(end,2) -q(end,1) q(end,4) -q(end,3);
q(end,1) q(end,2) q(end,3) q(end,4)];
E=[zeros(3) zeros(3) zeros(3,4) zeros(3,1);
zeros(3) eye(3) zeros(3,4) w(end,:)’;
zeros(4,3) zeros(4,3) Q_n’*Q_n zeros(4,1);
zeros(1,3) w(end,: ) zeros(1,4) w(end,:)*w(end,:)’];

for k=1:n-1
disp(k);
X(k,:)= [u(k,:), w(k+1,:), q(k+1,:)];%dim=nx10

R{k}=0.5*dt*[-q(k,4) q(k,3) -q(k,2);
    -q(k,3) -q(k,4) q(k,1);
    q(k,2) -q(k,1) -q(k,4);
    q(k,1) q(k,2) q(k,3)];

F{k}=[-dt*eye(3,3) J zeros(3,4);
    zeros(4,3) R{k} eye(4)];

y{k}=[J(1,1)*w(k,1)+dt*(J(2,2)-J(3,3))*w(k,2)*w(k,3);
    J(2,2)*w(k,2)+dt*(J(3,3)-J(1,1))*w(k,3)*w(k,1);
    J(3,3)*w(k,3)+dt*(J(1,1)-J(2,2))*w(k,1)*w(k,2);
    q(k,:)'];

cvx_begin
cvx_solver sdpt3
% cvx_precision medium

variables u(n-1,3) q(n,4) w(n,3)

% expression X(n-1,10)
% X(k,:)= [u(k,:), w(k+1,:), q(k+1,:)];%dim=nx10

minimize(sum_square_abs(w(k+1,:)'-w(end,:)')+sum_square_abs(Q_n*q(k+1,:)'))
%minimize(norm(w(k+1,:)'-w(end,:)')+norm(C*(Q_n*q(k+1,:)')))
%minimize ([X(k,:)';1]'*E*[X(k,:)';1]+norm(q(k+1,:)))
%minimize (quad_form([X(k,:)';1],E)+sum_square_abs(q(k+1,:)))

subject to

F{k}*X(k,:)' == y{k};

abs(G1*X(k,:)')<= gamma1*ones(3,1);
abs(G2*X(k,:)')<= gamma2*ones(3,1);

X(k,:)*H*X(k,:)' <=0;

cvx_end
Xstore(k,:)=[u(k,:) w(k+1,:) q(k+1,:)];

end


(Mark L. Stone) #2

Your program is not reproducible, as it is missing inputs and/or functions.

I will presume your results pertain to the version as provided, i.e. code in your CVX program which is commented out really is commented out.

It appears that in reality, contrary to your intention, each of the n-1 instances of invoking CVX are unconstrained minimization, for which the k+1 'th columns of optimal w and q are zeros, achieving optimal objective of 0. Your supposed constraints, which all involve X, are not actually CVX constraints, because they don’t involve any CVX variables or expressions; and they are just evaluated as logical (true or false) expressions by MATLAB, without involving or affecting CVX at all. And the CVX variable u is never used, so its value upon exit from CVX can be whatever CVX feels like.

If X is defined outside CVX, but not declared as a variable or expression or created as an expression inside CVX, then it is treated by CVX as just being input data (MATLAB variable). Anything already defined outside CVX, but declared as a variable or expression or created as an expression inside CVX will overwrite (supersede) any value it had prior to CVX.


#3

Thanks so much for an expedite respond. So, is it ok define variables as u(n-1,3) q(n,4) w(n,3) and X(n-1,10)? Also is it necessary to define the values of R, F, y, which are variable matrices as expressions? Additionally, can the constant matrices like G1, G2, H used for framing constraints be expressions or normal matlab matrices?


#4

Also is it ok to define:
variables u(n-1,3) q(n,4) w(n,3) X(n-1,10)
expression X(n-1,10)
for k=1:n-1
X(k,:)= [u(k,:), w(k+1,:), q(k+1,:)];%dim=nx10
end


(Mark L. Stone) #5

If each instance of CVX only involves a single value of k, then you are probably better off defining all CVX variables and expressions so as to eliminate that index, i.e., make u, q, w, and X be 1-D vectors. If you want, you should save all the optimal values into a MATLAB 2-D array (matrix), and update that after CVX completion each time through the loop.

If R, F, y are all created based on initial values (1st time through the for loop), or optimal; values from the previous time through the loop, then don’t declare them as CVX variables or expressions, because that would result in CVX over-writing their values, which I don’t think is what you want.

If something depends on or is a variable being optimized in that run of CVX, then it must be a CVX variable or expression (not all CVX expressions have to be declared). If it is input data to the optimization problem for that run of CVX, then it should have a numerical value prior to use, and should not be a CVX variable or expression.

I think you can benefit from re-reading the CVX Users’ Guide, with some of these things in mind.


#6

The code I provided is the complete code except u(1,:)= (J*wdot(1,:)’+[0 -w(1,3) w(1,2);w(1,3) 0 -w(1,1);-w(1,2) w(1,1) 0]Jw(1,:)’)’;%zeros(1,3);% calculated from initial values of q and w

where crossm(w(1,:)’=[0, -w(1,3), w(1,2);w(1,3), 0, -w(1,1);-w(1,2), w(1,1), 0]


#7

Dear Dr. Mark, Thanks so much for these valuable inputs. I will read it again religiously.


#8

Dear Dr. Mark, I modified my code keeping in mind the info you provided and user-manual’s info. Still, my objective value returns infinity for the k=1’st time of invoking CVX. I am not able to figure out why?
Can you please help. Here is the code.

clc;clear all;close all;

ti=0;tf=20;dt=1;t=[0:dt:tf]’;
n =((tf-ti)/dt);%divide t0-tf into n intervals

J=diag([100,200,300]);

W=zeros(n,3);
w=W(2:n,:);

Q=zeros(n,4);
q=Q(2:n,:);

% Initial and final conditions of angular velocity and quaternions
Q(1,:)=[0.5, 0.5, 0.5, 0.5];
W(1,:)=[0.0,0.0,0.0]*pi/180;

Q(end,:)=[0.0258, 0.0258, 0.999, 0.0258];
W(end,:)=[0,0,0];

G1=[eye(3,3) zeros(3,7)];%dim=3x10
gamma1=01;% control bound (rad/s^2)

G2=[zeros(3,3) eye(3) zeros(3,4)];%dim=3x10
gamma2=.3;% angular velocity bound (rad/s)

% Non-convex FoV exclusion constraint
xs_I=[0; 0; 1];
yc_B=[0.750;0.4330;0.500];
A =xs_I* yc_B’+yc_B * xs_I’-(xs_I’ * yc_B+cos(50 * pi/180)) * eye(3);
b = [0 -xs_I(3) xs_I(2);xs_I(3) 0 -xs_I(1);-xs_I(2) xs_I(1) 0] * yc_B;
d = xs_I’ * yc_B-cos(50 * pi/180);
Atilda = [A b;
b’ d];

eig_atilda=eig(Atilda);
tol=1e-12;
mu=max(tol,abs(min(eig(Atilda))));
Atilda1=Atilda+mu*eye(4);
eig_atilda1=eig(Atilda1);

H = [zeros(6,6) zeros(6,4);zeros(4,6) Atilda1];
C=[eye(3) zeros(3,1)];
G3=[zeros(4,3) zeros(4,3) eye(4,4)];
y=zeros(7,1,n-1);
R=zeros(4,3,n-1);
F=zeros(7,10,n-1);

Q_n=[Q(end, 4) Q(end,3) -Q(end,2) -Q(end,1);
-Q(end,3) Q(end,4) Q(end,1) -Q(end,2);
Q(end,2) -Q(end,1) Q(end,4) -Q(end,3);
Q(end,1) Q(end,2) Q(end,3) Q(end,4)];

for k=1:n-1
disp(k);

R(:,:,k)=0.5*dt*[-Q(k,4) Q(k,3) -Q(k,2);
    -Q(k,3) -Q(k,4) Q(k,1);
    Q(k,2) -Q(k,1) -Q(k,4);
    Q(k,1) Q(k,2) Q(k,3)];

F(:,:,k)=[-dt*eye(3,3) J zeros(3,4);
    zeros(4,3) R(:,:,k) eye(4)];

y(:,:,k)=[J(1,1)*W(k,1)+dt*(J(2,2)-J(3,3))*W(k,2)*W(k,3);
    J(2,2)*W(k,2)+dt*(J(3,3)-J(1,1))*W(k,3)*W(k,1);
    J(3,3)*W(k,3)+dt*(J(1,1)-J(2,2))*W(k,1)*W(k,2);
    Q(k,:)'];

cvx_begin

variables X(n-1,10)

minimize (sum_square_abs(G2*X(k,:)'-W(end,:)')+sum_square_abs(C*(Q_n*(G3*X(k,:)'))))
subject to

F(:,:,k)*X(k,:)' == y(:,:,k);
abs(G1*X(k,:)')<= gamma1*ones(3,1);
abs(G2*X(k,:)')<= gamma2*ones(3,1);
X(k,:)*H*X(k,:)' <=0;

cvx_end

Xstore(k,:)=X(k,:);
Q(k+1,:)=G3*X(k,:)';
W(k+1,:)=G2*X(k,:)';

end


(Mark L. Stone) #9

The problem provided to CVX when k = 1 is infeasible, hence cvx_optval = +Inf. I leave to you the burden of figuring out what your code and data should be for what you are trying to do.


#10

Hi Dr. Mark, is it valid to describe cell arrays in expression? For eg . expression W1{n-1};
for k=1:n-1
W1{k}=[-J(1,1) , C_w3(3*k,: )X, C_w3(3k-1,: )X;
C_w3(3
k,: )X, -J(2,2), C_w3(3k-2,: )X;
C_w3(3
k-1,: )X, C_w3(3k-2,: )*X , -J(3,3)];
end


(Mark L. Stone) #11

i don’t believe you can do that. You can’t declare cells as CVX variables or expressions, although you can put CVX variables and expressions in a cell.

I recommend you consider using multidimensional arrays, with k being the last of the 3 indices. So for a 3D array W1, W1(:,:,k) would be a 2-D matrix. See mcg 's answers at How to change the number of variables based on the iteration number? and Cvx sdp mode and cell arrays .


#12

Dear Dr. Mark, Thanks so much for all the insight into the CVX working. It has helped me a LOT as now my code is working and I understood CVX a little better while troubleshooting my code. Thanks again so much!