Help! Status: Unbounded

AandZERO · April 11, 2024, 8:32am

the following is my code, which is a SDP problem, however, the result shows “Status: Unbounded”.
I want to know why this happens.

clc
clear all
close all

%% 测量点位置
Radius=1.500;
Angle=linspace(0,2pi,30);
PosLoc=[Radiuscos(Angle(1:end-1));Radiussin(Angle(1:end-1))];
PosLoc=PosLoc’;
xShip=PosLoc(:,1)+1000.001randn(size(PosLoc(:,1)));
yShip=PosLoc(:,2)+1000.001*randn(size(PosLoc(:,1)));
zShip=0; % 发射换能器的入水深度

%% 目标位置
Range=600;[200:100:1000];
xTar=0.001Range/sqrt(2);
yTar=0.001Range/sqrt(2);
NN=length(Range); % 阵元数目
zTar=0;

%% 添加测量误差
detat=1*0.001; % 时延测量误差

c0=1.500;
detac=010.001; % 声速测量误差

detaa=010.001; % 节点位置的随机误差

T0=5;

%% 主程序
for nn=1:NN
nn

for mc=1:1
    
    %%%% 构造带误差的测量值
    xShipMea=xShip+detaa*randn(size(xShip));
    yShipMea=yShip+detaa*randn(size(yShip));
    
    tMea=sqrt((xShip-xTar(nn)).^2+(yShip-yTar(nn)).^2+(zShip-zTar).^2)/c0+detat*randn(size(xShip))+(0:length(xShip)-1)'*T0;
    
    c=c0+detac*randn;
    
    %%% 主程序
    s=[xShipMea(2:end)-xShipMea(1),yShipMea(2:end)-yShipMea(1)];
    
    t=tMea(2:end)-tMea(1);
    
    M=length(t);
    
    h1=sum(s.^2,2)-t.^2*c^2;
    
    A1=2*s;
    
    N=(1:M)';
    B1=[2*t*c,-2*N*c,-2*N.*t*c^2,N.^2*c^2];
    
    Qt=detat^2*(ones(length(tMea)-1)+eye(length(tMea)-1));        
    
    W=eye(M);
    for it=1
        
        F=[A1 B1];
        
        PHI=[F'*W*F,-F'*W*h1;
            -h1'*W*F,h1'*W*h1];
        
        cvx_begin sdp
        variable y(6,1);
        variable Y(6,6);
        minimize(trace([Y,y;y',1]*PHI));
        subject to
        
        [Y,y;y',1]==semidefinite(7);
        
        Y(1,1)+Y(2,2)-Y(3,3)==0;
        
        Y(3,5)-y(4)==0;
        
        Y(5,5)-y(6)==0;
        
        Y(3,6)-Y(4,5)==0;

% norm(Y(1,1)+Y(2,2)) <= y(3);
%
% norm(Y(1,5)+Y(2,5)) <= y(4);
%
% norm(Y(1,6)+Y(2,6))<=Y(3,6);

        cvx_end

        xFinal=y(1);
        yFinal=y(2);
        TFinal=y(5);
        
        C1=-2*diag(sqrt((xFinal-s(:,1)).^2+(yFinal-s(:,2)).^2));            
        W=eye(M)/(C1*Qt*C1');
    end
    
    xFinal=xFinal+xShipMea(1);
    yFinal=yFinal+yShipMea(1);
    
    xSDR(mc,nn)=xFinal;
    ySDR(mc,nn)=yFinal;
    
end

end

%%% 均方根误差
RmseSDR=1000*sqrt(sum((xSDR-ones(mc,1)*xTar).^2+(ySDR-ones(mc,1)*yTar).^2,1)/mc);

%%%% 估计偏差
BiasSDR=1000*sqrt((mean(xSDR-ones(mc,1)*xTar,1)).^2+(mean(ySDR-ones(mc,1)*yTar,1)).^2);

%% 绘图和保存数据
figure;
hold on
scatter(1000xShip,1000yShip,‘blue.’);
scatter(1000xTar,1000yTar,‘red.’);
legend(‘水面船航路’,‘海底阵阵元’)
xlabel(‘x(m)’);
ylabel(‘y(m)’);
xlim([-2000 2000])
ylim([-2000 2000])
axis equal
grid on
box on
name=‘几何态势图’;
% SavePicture( gcf,name,0 );

figure;
hold on
plot(Range,RmseSDR,‘red-*’,‘LineWidth’,1.5)
box on
grid on
set(gca,‘GridLineStyle’,‘–’)
xlabel(‘海底阵阵元和原点之间的水平距离(m)’)
ylabel(‘海底阵阵元位置的均方根误差(m)’)
name=‘均方根误差’;

figure;
hold on
plot(Range,BiasSDR,‘red-*’,‘LineWidth’,1.5)
box on
grid on
set(gca,‘GridLineStyle’,‘–’)
xlabel(‘海底阵阵元和原点之间的水平距离(m)’)
ylabel(‘海底阵阵元位置的估计偏差(m)’)
name=‘估计偏差’;
% SavePicture( gcf,name,0 );

Mark_L_Stone · April 11, 2024, 11:18am

Follow the advice in Debugging unbounded models - YALMIP

AandZERO · April 11, 2024, 1:12pm

I get rid of the object, and it shows “Solved”.
the object is reserved, and I add the constraint “trace([Y,y;y’,1]) >= 20”, it shows “Unbpunded”, the result “y” do not satisfy the constraint “trace([Y,y;y’,1]) >= 20”.

is there some problem with the object?

Mark_L_Stone · April 11, 2024, 1:37pm

I am surprised there is not a warning message for the semidefinite constraint due to not declaring Y symmetric. I think you should declare it symmetric, even though that might not be contributing to your difficulties.

Variable values are meaningless unless the problem is reported solved to optimality.

The problem will never be reported unbounded if there is no objective.

Perhaps you should concentrate on putting explicit bounds on all the variables while keeping the objective, and see which variables end up on those bounds. For example: -1000 <= Y(:) <= 1000, -1000 <= y <= 1000.

AandZERO · April 11, 2024, 2:06pm

I add some explicit bounds as follows

        for ii=3:6
            for jj=3:6
            0<=Y(ii,jj)<=1000
            end
            0<=y(ii)<=1000
        end

the result is “Solved”, but the answer is not optimal, I can find an input which satisfies the constraints and results in a smaller Optimal value.

Mark_L_Stone · April 11, 2024, 2:12pm

Show us the solver and CVX output, and show us CVX’s solution and your evidence that you have a better solution.

The purpose of adding the constraints is to help you diagnose why the problem is unbounded, and help you figure out what to do about it. If the problem is unbounded, at least one of the variables (elements) should be on a bound when bounds are added.

AandZERO · April 11, 2024, 2:21pm

the following is code:
clc
clear all
close all

%% 测量点位置
Radius=1.500;
Angle=linspace(0,2pi,60);
PosLoc=[Radiuscos(Angle(1:end-1));Radiussin(Angle(1:end-1))];
PosLoc=PosLoc’;
xShip=PosLoc(:,1)+10000.001randn(size(PosLoc(:,1)));
yShip=PosLoc(:,2)+10000.001*randn(size(PosLoc(:,1)));
zShip=0; % 发射换能器的入水深度

%% 目标位置
Range=600;[200:100:1000];
xTar=0.001Range/sqrt(2);
yTar=0.001Range/sqrt(2);
NN=length(Range); % 阵元数目
zTar=0;

%% 添加测量误差
detat=1*0.001; % 时延测量误差

c0=1.500;
detac=010.001; % 声速测量误差

detaa=010.001; % 节点位置的随机误差

T0=5;

%% 主程序
for nn=1:NN
nn

for mc=1
    
    %%%% 构造带误差的测量值
    xShipMea=xShip+detaa*randn(size(xShip));
    yShipMea=yShip+detaa*randn(size(yShip));
    
    tMea=sqrt((xShip-xTar(nn)).^2+(yShip-yTar(nn)).^2+(zShip-zTar).^2)/c0+0*detat*randn(size(xShip))+(0:length(xShip)-1)'*T0;
    
    c=c0+detac*randn;
    
    %%% 主程序
    s=[xShipMea(2:end)-xShipMea(1),yShipMea(2:end)-yShipMea(1)];
    
    t=tMea(2:end)-tMea(1);
    
    M=length(t);
    
    h1=sum(s.^2,2)-t.^2*c^2;
    
    A1=2*s;
    
    N=(1:M)';
    B1=[2*t*c,-2*N*c,-2*N.*t*c^2,N.^2*c^2];
    
    Qt=detat^2*(ones(length(tMea)-1)+eye(length(tMea)-1));        
    
    W=eye(M);
    for it=1
        
        F=[A1 B1];
        
        PHI=[F'*W*F,-F'*W*h1;
            -h1'*W*F,h1'*W*h1];
        
        cvx_begin sdp
        variable y(6,1);
        variable Y(6,6) symmetric;
        minimize(trace([Y,y;y',1]*PHI));
        subject to
        
        [Y,y;y',1]==semidefinite(7);
        
        Y(1,1)+Y(2,2)-Y(3,3)==0;
        
        Y(3,5)-y(4)==0;
        
        Y(5,5)-y(6)==0;
        
        Y(3,6)-Y(4,5)==0;
        
        for ii=3:6
            for jj=3:6
            0<=Y(ii,jj)<=1000
            end
            0<=y(ii)<=1000
        end

% norm(Y(1,1)+Y(2,2)) <= y(3);
%
% norm(Y(1,5)+Y(2,5)) <= y(4);
%
% norm(Y(1,6)+Y(2,6))<=Y(3,6);

        cvx_end

        xFinal=y(1);
        yFinal=y(2);
        TFinal=y(5);

% C1=-2diag(sqrt((xFinal-s(:,1)).^2+(yFinal-s(:,2)).^2));
% W=eye(M)/(C1Qt*C1’);
end

    xFinal=xFinal+xShipMea(1);
    yFinal=yFinal+yShipMea(1);
    
    xSDR(mc,nn)=xFinal;
    ySDR(mc,nn)=yFinal;
    
    
    ux=xTar(1)-xShipMea(1);
    uy=yTar(1)-yShipMea(1);
    
    y0=[ux;uy;sqrt(ux^2+uy^2);sqrt(ux^2+uy^2)*T0;T0;T0^2];
    Y0=y0*y0';
    trace([Y0,y0;y0',1]*PHI)
    
end

end

%%% 均方根误差
RmseSDR=1000*sqrt(sum((xSDR-ones(mc,1)*xTar).^2+(ySDR-ones(mc,1)*yTar).^2,1)/mc);

%%%% 估计偏差
BiasSDR=1000*sqrt((mean(xSDR-ones(mc,1)*xTar,1)).^2+(mean(ySDR-ones(mc,1)*yTar,1)).^2);

%% 绘图和保存数据
figure;
hold on
scatter(1000xShip,1000yShip,‘blue.’);
scatter(1000xTar,1000yTar,‘red.’);
legend(‘水面船航路’,‘海底阵阵元’)
xlabel(‘x(m)’);
ylabel(‘y(m)’);
xlim([-2000 2000])
ylim([-2000 2000])
axis equal
grid on
box on
name=‘几何态势图’;
% SavePicture( gcf,name,0 );

figure;
hold on
plot(Range,RmseSDR,‘red-*’,‘LineWidth’,1.5)
box on
grid on
set(gca,‘GridLineStyle’,‘–’)
xlabel(‘海底阵阵元和原点之间的水平距离(m)’)
ylabel(‘海底阵阵元位置的均方根误差(m)’)
name=‘均方根误差’;

figure;
hold on
plot(Range,BiasSDR,‘red-*’,‘LineWidth’,1.5)
box on
grid on
set(gca,‘GridLineStyle’,‘–’)
xlabel(‘海底阵阵元和原点之间的水平距离(m)’)
ylabel(‘海底阵阵元位置的估计偏差(m)’)
name=‘估计偏差’;
% SavePicture( gcf,name,0 );

the following is result, the optimal value of CVX is 1338.42
y0 and Y0 is a better input and its object value 1.2730e-04

Mark_L_Stone · April 11, 2024, 2:28pm

Perhaps you can explicitly and directly assign numerical values to y0 and Y0, then show us the constraints are satisfied.

In any event, given the huge objective value as experienced by SDPT3 (4.314e11), probably all results are suspect. Perhaps bad numerical scaling. Perhaps that is a consequence of, or at least exacerbated by, the large bounds put on the problem. You need to figure out what additional constraints are needed, or perhaps that the value of PHI is not compatible with the problem being bounded.

AandZERO · April 12, 2024, 8:43am

for uu=1
ux=xTar(1)-xShipMea(1);
uy=yTar(1)-yShipMea(1);
ru=sqrt(ux^2+uy^2);
y0=[ux;uy;ru;ru**T0;T0;T0^2];
Y0=y0*y0’;
end

the above is the more accurate input (y0,Y0) whose object value is 1.2730e-04, and it satifies the following constraint

[Y,y;y’,1]==semidefinite(7);

    Y(1,1)+Y(2,2)-Y(3,3)==0;

where Y0(1,1)+Y0(2,2)-Y0(3,3)=y0(1)*y0(1)+y0(2)*y0(2)-y0(3)*y0(3)=ux^2+uy^2-ru^2=0

    Y(3,5)-y(4)==0;    where Y0(3,5)=y0(3)*y0(5)=ru*T0=y0(4)
    
    Y(5,5)-y(6)==0;     where Y0(5,5)=y0(5)*y0(5)=T0*T0=y0(6)
    
    Y(3,6)-Y(4,5)==0;   where Y0(3,6)=y0(3)*y0(6)=ru*T0^2=(ru*T0)*T0=y0(4)*y0(5)=Y0(4,5)

Mark_L_Stone · April 12, 2024, 10:59am

You haven’t provided Y0 and y0 in an explicit and direct form we can reproduce. For instance,. y0=[.18534 .232233 33333.2 3.45566 1.234556 .343434343];

In any event which elements of the "optimized’ Y0, t0 are on bounds? What is the deficiency of your model which causes that? Perhaps PHI is causing that?

AandZERO · April 12, 2024, 11:53am

Every time the program runs, the values of the variables are different, and I cannot give a fixed set of numbers. the following is my matlab code , the variables y0 and y0 can be found

clc
clear all
close all

format long

for zzz=1

%% 测量点位置
Radius=1.500;
Angle=linspace(0,2*pi,60);
PosLoc=[Radius*cos(Angle(1:end-1));Radius*sin(Angle(1:end-1))];
PosLoc=PosLoc';
xShip=PosLoc(:,1)+1000*0.001*randn(size(PosLoc(:,1)));
yShip=PosLoc(:,2)+1000*0.001*randn(size(PosLoc(:,1)));
zShip=0;  % 发射换能器的入水深度

%% 目标位置
Range=600;[200:100:1000];
xTar=0.001*Range/sqrt(2);
yTar=0.001*Range/sqrt(2);
NN=length(Range);  % 阵元数目
zTar=0;

%% 添加测量误差
detat=1*0.001;  % 时延测量误差

c0=1.500;
detac=0*1*0.001;  % 声速测量误差

detaa=0*1*0.001;  %  节点位置的随机误差

T0=5;

%% 主程序
for nn=1:NN
    nn
    
    for mc=1
        
        %%%% 构造带误差的测量值
        xShipMea=xShip+detaa*randn(size(xShip));
        yShipMea=yShip+detaa*randn(size(yShip));
        
        tMea=sqrt((xShip-xTar(nn)).^2+(yShip-yTar(nn)).^2+(zShip-zTar).^2)/c0+0*detat*randn(size(xShip))+(0:length(xShip)-1)'*T0;
        
        c=c0+detac*randn;
        
        %%% 主程序
        s=[xShipMea(2:end)-xShipMea(1),yShipMea(2:end)-yShipMea(1)];
        
        t=tMea(2:end)-tMea(1);
        
        M=length(t);
        
        h1=sum(s.^2,2)-t.^2*c^2;
        
        A1=2*s;
        
        N=(1:M)';
        B1=[2*t*c,-2*N*c,-2*N.*t*c^2,N.^2*c^2];
        
        Qt=detat^2*(ones(length(tMea)-1)+eye(length(tMea)-1));
        
        W=eye(M);
        
        ux=xTar(1)-xShipMea(1);
        uy=yTar(1)-yShipMea(1);
        
        for it=1
            
            F=[A1 B1];
            
            PHI=(([F -h1])'*W*([F -h1]));
            
            cvx_begin sdp
            variable y(6,1);
            variable Y(6,6) symmetric;
            minimize(trace([Y,y;y',1]*PHI));
            subject to
            
            [Y,y;y',1]==semidefinite(7);
            
            Y(1,1)+Y(2,2)-Y(3,3)==0;
            
            Y(3,5)-y(4)==0;
            
            Y(5,5)-y(6)==0;
            
            Y(3,6)-Y(4,5)==0;
            
            for ii=3:6
                for jj=3:6
                    0<=Y(ii,jj)<=1000;
                end
            end
            
            
            cvx_end
            
            xFinal=y(1);
            yFinal=y(2);
            TFinal=y(5);
            
        end
        
        xFinal=xFinal+xShipMea(1);
        yFinal=yFinal+yShipMea(1);
        
        xSDR(mc,nn)=xFinal;
        ySDR(mc,nn)=yFinal;
        
        y0=[ux;uy;sqrt(ux^2+uy^2);sqrt(ux^2+uy^2)*T0;T0;T0^2];
        Y0=y0*y0';
        trace([y0*y0',y0;y0',1]*PHI)
        
    end
    
    
end

end

AandZERO · April 12, 2024, 11:57am

up to now, i attempt to add some constraints on y0 and Y0, and i do not find some variable on bounds

for example
for ii=3:6
for jj=3:6
0<=Y(ii,jj)<=1000;
end
end

Mark_L_Stone · April 12, 2024, 12:11pm

You are generating random input data using randn. You could pick one value of those inputs. Then the problem should be reproducible and provide the same output whenever it is run.

It appears to me that PHI is the cause of your difficulties. Maybe you should concentrate on that. For instance, if PHI were the Identity matrix, I don’t believe the problem would be unbounded.