An iterative alternating optimization algorithm on Mosek

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

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, = 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

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.

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.

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.

Thank you! I’ll learn it carefully!

Hi, have you solved this problem? I come up with a similar issue.