An iterative alternating optimization algorithm on Mosek

1

I have an iterative alternating optimization algorithm on Mosek. It runs successfully in initial several times, but will be the Infeasible after. Strangly, the iteration number before infeasible status increases, when I increase a parameter (which is positive correlation with the size of variable) in my environment.
I don’t know why it will be infeasible several times after. Can you help me ?
The print of solver is：

Calling Mosek 9.1.9: 398 variables, 172 equality constraints
For improved efficiency, Mosek is solving the dual problem.

MOSEK Version 9.1.9 (Build date: 2019-11-21 11:34:40)
Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com
Platform: Windows/64-X86

Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 172
Cones : 72
Scalar variables : 398
Matrix variables : 0
Integer variables : 0

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 36
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries : 2 time : 0.00
Lin. dep. - tries : 1 time : 0.01
Lin. dep. - number : 0
Presolve terminated. Time: 0.03
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 172
Cones : 72
Scalar variables : 398
Matrix variables : 0
Integer variables : 0

Optimizer - threads : 6
Optimizer - solved problem : the primal
Optimizer - Constraints : 108
Optimizer - Cones : 72
Optimizer - Scalar variables : 323 conic : 216
Optimizer - Semi-definite variables: 0 scalarized : 0
Factor - setup time : 0.00 dense det. time : 0.00
Factor - ML order time : 0.00 GP order time : 0.00
Factor - nonzeros before factor : 2736 after factor : 2736
Factor - dense dim. : 0 flops : 3.13e+05
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 1.1e+00 2.4e+03 2.2e+00 0.00e+00 1.248000000e+00 0.000000000e+00 1.0e+00 0.05
1 4.9e-01 1.1e+03 1.5e+00 -9.97e-01 -7.025453946e+00 -7.052587232e+00 4.5e-01 0.13
2 2.4e-01 5.3e+02 1.0e+00 -9.89e-01 -1.840778278e+01 -1.626236371e+01 2.2e-01 0.13
3 8.2e-02 1.8e+02 5.7e-01 -9.60e-01 -6.758053191e+01 -5.805701666e+01 7.4e-02 0.13
4 7.9e-03 1.7e+01 1.1e-01 -8.10e-01 -5.241792360e+02 -4.771390819e+02 7.2e-03 0.14
5 1.8e-03 3.9e+00 2.3e-02 6.37e-03 -7.653769814e+02 -7.276327213e+02 1.6e-03 0.14
6 4.5e-04 9.9e-01 1.2e-02 -3.15e-01 -8.729409703e+02 -7.086695354e+02 4.1e-04 0.14
7 3.6e-06 7.8e-03 8.7e-04 -8.11e-01 -1.482981838e+04 -7.537200910e+02 3.3e-06 0.14
8 8.8e-14 1.7e-10 2.2e-07 -9.99e-01 -9.024137793e+11 -7.015536776e+02 7.1e-14 0.14
Optimizer terminated. Time: 0.19

Interior-point solution summary
Problem status : DUAL_INFEASIBLE
Solution status : DUAL_INFEASIBLE_CER
Primal. obj: -2.3809033638e-01 nrm: 2e+05 Viol. con: 7e-13 var: 0e+00 cones: 0e+00
Optimizer summary
Optimizer - time: 0.19
Interior-point - iterations : 8 time: 0.16
Basis identification - time: 0.00
Primal - iterations : 0 time: 0.00
Dual - iterations : 0 time: 0.00
Clean primal - iterations : 0 time: 0.00
Clean dual - iterations : 0 time: 0.00
Simplex - time: 0.00
Primal simplex - iterations : 0 time: 0.00
Dual simplex - iterations : 0 time: 0.00
Mixed integer - relaxations: 0 time: 0.00

Status: Infeasible
Optimal value (cvx_optval): +Inf

2

LEVEL=2;EPSILON=10;V=3;K=25;
M=1e6;
Q_INI =50;Q_NUM=50;Q_STEP = 5;
q_lst = (Q_NUM-Q_INI)/Q_STEP+1;
B = 20e7;
P_THR = 1e-23; C = 0;
alpha_int = 250;
weight=5;
ALPHA_P=100000;
alpha_cm_a=15e-7weight;alpha_cp_a=4;
alpha_cm_g=10e-7weight;alpha_cp_g=4;
cn =[150, 150, 150, 150, 150, 150, 150]*3;% cn = round(rand(1,7)30)+100;
pmax = [0.1,0.1,0.1]; PMAX=0.1; cigma = 3.981e-21;
uav_num = 3;bs_num = 4;node_num = 7;
lg = zeros(bs_num,7)+1.5e8; lg=lg./M;
% lg2 = zeros(bs_num,uav_num)+1.5e8./2;
% LG1 =1.5e8; LG2 = 1.5e8;
tq = [10,13,10,20];
ser = [3,3,3;1,1,1;1,1,1;3,3,3].*5;data = [8,4;20,10;6,3;24,12].*25e6;rq = [1500,2500,1000,2500];
uav_pos = [100,300,50; 300,100,50;500,300,50]; % 5 6 7
bs_pos = [200,200,0; 200,400,0;% |1-3|
400,200,0;400,400,0];% |2-4|
node_pos = [bs_pos;uav_pos];
H = zeros(uav_num,node_num);
for i = 1:uav_num
for j = 1:node_num
H(i,j) = norm(uav_pos(i,:)-node_pos(j,:));
if H(i,j)~=0
H(i,j)=1/H(i,j)^3.72;
end
end
end
for qnum = 1:q_lst
q_num = Q_INI+(qnum-1)*Q_STEP;
ld = zeros(1,3);
iter_r=zeros(1,K);
for level = 1:LEVEL
q_list =1:q_lst;
[sd,q]=service_generate(q_num,1);
rqi = zeros(1,q_num);
cqi = zeros(q_num,3);
tqi = zeros(q_num,1);
lq = zeros(q_num,2);
delta_tq = zeros(q_num,3);
lkk = zeros(node_num, node_num, Q_NUM);
for qi = 1:q_num
type_qi = q(qi);
rqi(qi)=rq(type_qi);
cqi(qi,:)=ser(type_qi,:);
tqi(qi,1)=tq(type_qi);
lq(qi,:slight_smile: = data(type_qi,:);
end
link0=lq./tqi.*2;
link0=link0./M;
link1=link0; link2=link0;
lq=lq./M;
%% part-1
cvx_clear
zvq = ones(q_num,1);r_k = zeros(q_lst,1);
z_k = zeros(q_num,K);p_k = zeros(uav_num, node_num,K);
b_k = zeros(uav_num, node_num,K);c0_k = zeros(node_num,K);
this_link_k=zeros(node_num,node_num,K);
r_k2 = zeros(q_lst,1);
z_k2 = zeros(q_num,K);p_k2 = zeros(uav_num, node_num,K);
b_k2 = zeros(uav_num, node_num,K);c0_k2 = zeros(node_num,K);
this_link_k2=zeros(node_num,node_num,K);
this_target=zeros(q_num,K);
k=1;flag=0;
while (k<=K && flag==0)
%% CVX-3-1
cvx_begin
variable x(3, node_num, q_num) nonnegative;
variable y(q_num, node_num, node_num, 2) nonnegative;
variable vp(uav_num, node_num) nonnegative;
variable vb(uav_num, node_num) nonnegative;
variable z(q_num) nonnegative;
expressions c0(node_num);
expressions lad(3);
expressions p(uav_num, node_num) nonnegative;
expressions this_link(7,7);
expressions b(uav_num, node_num) nonnegative;
b=vb.*1e6;
p=vp./1000;
for qi = 1:q_num
if (link2(qi,1)==0 || link2(qi,2)==0 || sum(ismember(isnan(link2),1),‘all’)>=1)
link2(qi,:)=lq(qi,:)./tqi(qi,:).*2;
end
for n = 1:node_num
if qi==1
c0(n) = cqi(qi,:)*x(:,n,qi);
else
c0(n) = cqi(qi,:)*x(:,n,qi) +c0(n);
end
end
end
for n = 1:bs_num
for m = 1:node_num
this_link(n,m)=y(:,n,m,1)'*link2(:,1)+y(:,n,m,2)'link2(:,2);
end
end
for n = 1:uav_num
for m = 1:node_num
this_link(n+4,m)= y(:,n+4,m,1)'link2(:,1)+y(:,n+4,m,2)'link2(:,2);
end
lad(n)=alpha_cp_ac0(n+4)+ALPHA_Psum(p(n,:));
end
maximize sfc_reward(rqi,z,M.this_link,c0,p)-alpha_int(sum(z)-sum(zvq.^2)-2.zvq’(z-zvq))
subject to
y(:)<=1;
z(:)<=1;
for qi = 1:q_num
x(1, sd(qi,1), qi)==z(qi); % constrain-1
x(3, sd(qi,2), qi)==z(qi); % constrain-2
sum(y(qi,sd(qi,1),:,1))<=1;
sum(y(qi,:,sd(qi,2),2))<=1;
for n = 1:node_num
for fi = 1:2
sum(y(qi,n,:,fi))-sum(y(qi,:,n,fi)) == x(fi, n, qi)-x(fi+1, n, qi);
end
end
sum(x(:,:,qi),1)==z(qi); % constrain-5
sum(x(:,:,qi),2)>=z(qi);
end
for n =1:node_num
c0(n)<=cn(n); % constrain 7
end
sum(vp,2)./1000<= PMAX;% constrain-9(None 8)
for n = 1:bs_num
for m = 1:node_num
this_link(n,m)<=lg(n,m); % constrain-11-a
end
end
% constrain-12 xf(y/x) = -rel_entr(x,x+y)
for n = 1:uav_num
for m = 1:node_num
this_link(n+4,m)<=-rel_entr(b(n,m), b(n,m) + p(n,m)H(n,m)/cigma)/log(2)./M;
end
end
for i =1:node_num
y(:,i,i,:)<=0;
end
sum(sum(b))<=B; % constrain-13
% constrain-8
(max(lad)-min(lad))<=EPSILONsqrt(3/2);
cvx_end
z=abs(z);
r_k2(k)=sfc_reward(rqi,z,this_link.M,c0,p)-alpha_int(sum(z)-sum(zvq.^2)-2.zvq’(z-zvq));

        %% part-2
        cvx_begin
        variable l(q_num,2) nonnegative;
        expressions sumt(q_num);
        expressions this_link(7,7);
        t0=cvx(zeros(q_num,1));
        t1=cvx(zeros(q_num,1));
        for n = 1:bs_num
            for m = 1:node_num
                this_link(n,m)=y(:,n,m,1)'*l(:,1)+y(:,n,m,2)'*l(:,2);
            end
        end
        for n = 1:uav_num
            for m = 1:node_num
                this_link(n+4,m)=y(:,n+4,m,1)'*l(:,1)+y(:,n+4,m,2)'*l(:,2);
            end
        end
        for qi = 1:q_num
            if z(qi)<=0.1
                continue
            else
                t0(qi,1)=lq(qi,1)*inv_pos(l(qi,1));
                t1(qi,1)=lq(qi,2)*inv_pos(l(qi,2));                    
                sumt(qi)=t0(qi)+t1(qi);
            end
        end
        minimize (sum(sum(this_link(1:4,:)))+5*sum(sum(this_link(5:7,:))))
        subject to
        for qi = 1:q_num
            if z(qi)<=0.1
                
                continue
            else
                sumt(qi)<=tqi(qi);
            end
        end
        for n = 1:bs_num
            for m = 1:node_num
                this_link(n,m)<=lg(n,m);
            end
        end
        for n = 1:uav_num
            for m = 1:node_num
                this_link(n+4,m)<=-rel_entr(b(n,m), b(n,m) + p(n,m)*H(n,m)/cigma)/log(2)/M; % constrain-12   x*f(y/x) = -rel_entr(x,x+y)
            end
        end
        cvx_end
        if sum(ismember(isnan(l),1),'all')>=1
           1==1; 
        end
        link2=l;
        r_k(k)=sfc_reward(rqi,z,this_link.*M,c0,p)-alpha_int*(sum(z)-sum(zvq.^2)-2.*zvq'*(z-zvq));
        zvq = z;
        k=k+1;
    end % k
end

end

3

I would follow the advice at

to figure out why your problem is primal infeasible. It is primal infeasible because the dual problem is dual infeasible.

4

IN addition to the link posted by @Erling , perhaps the iterative alternating optimization algorithm as you have implemented it is unstable, and doesn’t necessarily converge to anything; and at some iteration, produces an infeasible problem. You can follow the advice in the link and carefully examine the input data to diagnose what is happening on that iteration.

5

Thank you for helping me. Actually, I find that PART-1 will be solved all the time and the INFEASIBLE usually appears in PART-2. And I types the result of PART-1 and find that it doesn’t subject to my constrains in both PART-1 and PART-2.
Additionaly, I find that the the expressions used in cvx must be definded before SUBJECT TO in cvx program and there may be absent in User Guide.

6

Thank you! I’ll learn it carefully!