I am trying to replicate the paper:https://asp-eurasipjournals.springeropen.com/track/pdf/10.1186/s13634-016-0326-2
In the paper it clearly defines the problem in equtation 45 and the constraint by equation 29, 30 and 44. But when I tried to implement I get the above error… My eintire code and the out put of the CVX are below:
Any comment please.
CODE:
N = 4; % filter order
Rs = 40;% stopband attenuation dB
Wc = 0.2; % centre of stopband
Wb = 1/128; % bandwidth of stopband. For 128x oversampling Wb = 1/128
Ws = Wc +0.5Wb[ -1, 1];% stopband edge frequency (normalised)
Wf = Wc + 0.3Wb[-1, 1]; % output filter passband
% [z, p,k ] = cheby2( N, Rs, Ws);
[z, p,k ] = butter( N, Ws);
[sos,g] = zp2sos(z,p,k);
[Heb, Hea ]=sos2tf(sos,g);
figure(656)
plot(abs(freqz(Heb,Hea,2001)))
hold on
%%%%%%
orderN=N;
filLen=2*orderN;
FirOrder=32;
jita=1.5;
%%%% constants
Ar=zeros(FirOrder,FirOrder);
for kt=1:FirOrder-1
Ar(kt,kt+1)=1;
end
Br=zeros(FirOrder,1);
Br(FirOrder)=1;
Dr=1;
%%%% design parameter
Ah=zeros(filLen,filLen);
for kt=1:filLen
Ah(filLen,kt)=-Hea(filLen-kt+2);
end
% Ah(filLen,filLen)=-1;
% Ah(filLen,1)=-1;
Bh=zeros(filLen,1);
Br(filLen)=1;
Ch=zeros(1,filLen);
for kt=1:filLen
Ch(kt)=Heb(filLen-kt+2)-Hea(filLen-kt+2)*Heb(1);
end
Dh=Heb(1);
% Ch(filLen)=1;
[Ah1,Bh1,Ch1,Dh1]=tf2ss(Heb,Hea);
%%%% derived parameter;
A=[Ar, BrCh;zeros(filLen,FirOrder),Ah];
B=[BrDh;Bh];
% Cr=zeros(1,filLen);
D=Dh;
%%%%Optimize to find out parameter
cvx_begin % quiet
variables P(filLen+FirOrder,FirOrder+filLen)
variables Pr(FirOrder,FirOrder)
variables Cr(1,FirOrder);
variables u2
P == semidefinite(FirOrder+filLen);
Pr == semidefinite(FirOrder);
%%%% Equation 29
M1=[P,PA,PB];
M2=[A’P,P,zeros(FirOrder+filLen,1)];
M3=[B’*P,zeros(1,FirOrder+filLen),1];
mof29=[M1;M2;M3];
%%%% Equation 30
C=[Cr,Ch];
M1=[u2,C,D];
M2=[C’,P,zeros(FirOrder+filLen,1)];
M3=[D’,zeros(1,FirOrder+filLen),1];
mof30=[M1;M2;M3];
%%% equation 44
M1=[-Pr,PrAr, Pr*Br,zeros(FirOrder,1)];
M2=[Ar’*Pr, -Pr,zeros(FirOrder,1),Cr’];
M3=[Br’*Pr,0,-jita^2,ones(1,FirOrder)];
M4=[zeros(1,FirOrder),Cr,1,-1];
mof44=[M1;M2;M3;M4];
minimize (u2);
subject to
mof29>=0;
mof30 >= 0;
mof44 <= 0;
% Cr(filLen)==1;
cvx_end
OUTPUT OF CVX
Calling SDPT3 4.0: 7247 variables, 1419 equality constraints
For improved efficiency, SDPT3 is solving the dual problem.
num. of constraints = 1419
dim. of sdp var = 72, num. of sdp blk = 2
dim. of linear var = 5860
dim. of free var = 39
2978 linear variables from unrestricted variable.
*** convert ublk to lblk
SDPT3: Infeasible path-following algorithms
version predcorr gam expon scale_data
HKM 1 0.000 1 0
it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime
0|0.000|0.000|5.2e+04|1.5e+03|3.1e+07| 1.722389e+02 0.000000e+00| 0:0:00| spchol 1 1
1|0.039|0.179|5.0e+04|1.3e+03|7.8e+07| 1.764777e+02 0.000000e+00| 0:0:00| chol 1 1
2|0.260|0.116|3.7e+04|1.1e+03|6.3e+07| 2.603584e+02 0.000000e+00| 0:0:00| chol 1 1
3|0.390|0.517|2.3e+04|5.4e+02|4.6e+07| 3.570574e+02 0.000000e+00| 0:0:01| chol 1 1
4|0.289|0.367|1.6e+04|3.4e+02|3.4e+07| 3.911117e+02 0.000000e+00| 0:0:01| chol 1 1
5|0.016|0.045|1.6e+04|3.3e+02|3.4e+07| 3.978931e+02 0.000000e+00| 0:0:01| chol 1 1
6|0.486|0.637|8.1e+03|1.2e+02|2.1e+07| 5.245779e+02 0.000000e+00| 0:0:01| chol 1 1
7|0.617|0.203|3.1e+03|9.5e+01|1.1e+07| 4.973036e+02 0.000000e+00| 0:0:01| chol 1 1
8|0.645|0.745|1.1e+03|2.4e+01|4.1e+06|-7.404704e+00 0.000000e+00| 0:0:01| chol 1 1
9|0.598|0.159|4.4e+02|2.2e+01|2.3e+06|-2.147230e+04 0.000000e+00| 0:0:02| chol 1 1
10|0.504|0.016|2.2e+02|2.5e+01|1.8e+06|-3.656464e+07 0.000000e+00| 0:0:02| chol 1 1
11|0.331|0.002|1.5e+02|2.9e+01|6.9e+08|-3.929227e+12 0.000000e+00| 0:0:02| chol 1 1
12|0.001|0.000|1.5e+02|3.2e+01|1.7e+12|-9.860969e+16 0.000000e+00| 0:0:02| chol 1 1
13|0.000|0.000|1.5e+02|3.6e+01|1.1e+15|-8.924573e+19 0.000000e+00| 0:0:02|
sqlp stop: primal or dual is diverging, 1.9e+16
number of iterations = 13
Total CPU time (secs) = 2.17
CPU time per iteration = 0.17
termination code = 3
DIMACS: 1.5e+02 0.0e+00 4.1e+01 0.0e+00 -1.0e+00 1.3e-05
Status: Failed
Optimal value (cvx_optval): NaN