Standard deviation can't be output.

clear;clc;close all;
%% Parameters
p=5;
s=5;
t=10;
for lambda=0.2:0.2:0.2
    for n=1:5
        pc=80+10*ones(p,s,t);%reshape(xlsread('data.xlsx','sheet1','A2:A33'),p,s,t);
        pc1=120+10*ones(p,s,t);
        oc=200+150*ones(p,s,t);%reshape(xlsread('data.xlsx','sheet1','B2:B33'),p,s,t);
        hc=5+5*ones(p,t);%reshape(xlsread('data.xlsx','sheet1','G2:G17'),p,t);
        tc=10+8*ones(s,t);%xlsread('data.xlsx','sheet1','D2:D5');
        D0=500+100*ones(p,t);%reshape(xlsread('data.xlsx','sheet1','H2:H17'),p,t);
        sc=800+400*ones(p,s,t);%reshape(xlsread('data.xlsx','sheet1','C2:C33'),p,s,t);
        eo=2+ones(p,s,t);%reshape(xlsread('data.xlsx','sheet1','P2:P33'),p,s,t);
        eh=1+2*ones(p,t);%reshape(xlsread('data.xlsx','sheet1','N2:N17'),p,t);
        ev=2+3*ones(p,s,t);%reshape(xlsread('data.xlsx','sheet1','I2:I33'),p,s,t);
        w=350+200*ones(1,t);%xlsread('data.xlsx','sheet1','O2:O5');
        dis=100+50*ones(1,s);%[120 150];
        cc=4000;
        r=[0.01,0.1,0.1,0.05,0.15];
        cp=50;%60+10*rand;
        b=2+0.3*rand; gam=2;
        eta1=500000; epsilon=0.9; eta=0.1;theta=0.1;%lambda=0.9;
        theta1=150;p0=[0.3;0.6;0.1];e0=[1;2;3];e=e0/norm(e0,1);
        %p0=xlsread('data.xlsx','box','ah2:ah81');e0=xlsread('data.xlsx','box','z2:z81');e=e0/norm(e0,1);
        m=3;
        x=sdpvar(p,s,t);
        y=binvar(p,s,t);
        x1=sdpvar(p,s,t);
        y1=binvar(p,s,t);
        z=sdpvar(p,t);
        sh=sdpvar(p,t);
        u=sdpvar(1);beta=sdpvar(1);beta0=sdpvar(1);tao=sdpvar(m,1);tao0=sdpvar(m,1);
        r1=sdpvar(1);r2=sdpvar(m,1);r3=sdpvar(m,1);r4=sdpvar(m,1);
        lambda1=sdpvar(1);lambda2=sdpvar(m,1);lambda3=sdpvar(m,1);lambda4=sdpvar(m,1);
        r10=sdpvar(1);r20=sdpvar(m,1);r30=sdpvar(m,1);r40=sdpvar(m,1);
        lambda10=sdpvar(1);lambda20=sdpvar(m,1);lambda30=sdpvar(m,1);lambda40=sdpvar(m,1);
        %% Objective function
        obj=lambda*beta+(lambda/(1-epsilon))*(tao'*p0+(lambda2'+lambda3')*e*gam)+(1-lambda)*(sum(sum(sum(eta1.*D0)))*ones(1,m)*p0...
            +(r2'+r3')*e*gam)+10*sum(squeeze(sum(sum((x+x1),1).*tc,3)).*dis) + sum(sum(sum(oc.*(y+y1))))+sum(sum(hc.*z))+u*cp...
            +sum(theta1.*r.*(sum(sum(y+y1,1),3)))+sum(sum(sum((x1).*pc1)))+sum(sum(eta1.*(-sum(x,2)-sum(x1,2))))+sum(z(:,10))...
            +sum(sum(sum((x1).*pc1)))+sum(sum(sum((x).*pc)));%
        
        %% Constraint
        C=[x(:)>=0];
        C=[C,x1(:)>=0];
        C=[C,z(:)>=0];
        C=[C,sh>=0];
        %% disruption
        for i=1:2
            y(:,i,:)=0;
        end
        C=[C,r2(:)>=0];C=[C,r3>=0];C=[C,r4(:)>=0];
        C=[C,lambda2(:)>=0];C=[C,lambda3>=0];C=[C,lambda4(:)>=0];
        C=[C,r20(:)>=0];C=[C,r30>=0];C=[C,r40(:)>=0];
        C=[C,lambda20(:)>=0];C=[C,lambda30>=0];C=[C,lambda40(:)>=0];
        C=[C,tao(:)>=0];C=[C,tao0(:)>=0];
        C=[C,x<=sc.*y];
        C=[C,x1<=(sc-x).*y1];
        C =[C,sum(z,1)<=w];
        for i=1:p
            for j=1:t
                C =[C,sum(r.*y1(i,:,j),1)<=0.2];
            end
        end
        C =[C,-squeeze(sum(x(:,:,1),2))-squeeze(sum(x1(:,:,1),2))+D0(:,1)+z(:,1)>=0];
        for i=2:t
            C =[C,-z(:,i-1) - squeeze(sum(x(:,:,i),2)) - squeeze(sum(x1(:,:,i),2))+D0(:,i)+z(:,i)>=0];
        end
        
        C=[C,sum(sum(sum(eta1.*D0)))==r1*e-r2+r3];
        C=[C,tao==lambda1*e-lambda2+lambda3];
        C=[C,tao>=sum(sum(sum(eta1.*D0)))-beta*e];
        
        for i=1:p
            for j=1:t
                C=[C,lambda*beta0+(lambda/(1-epsilon))*(tao0'*p0+(lambda20'+lambda30')*e*gam)+(1-lambda)*(theta*D0(i,j)*ones(1,m)*p0...
                    +(r20'+r30')*e*gam)>=eta*sum(x(i,:,j)+x1(i,:,j),2)];
            end
        end
        for i=1:p
            for j=1:t
                C=[C,theta*D0(i,j)==r10*e-r20+r30];
            end
        end
        C=[C,tao0==lambda10*e-lambda20+lambda30];
        for i=1:p
            for j=1:t
                C=[C,tao0>=theta*D0(i,j)-beta0*e];
            end
        end
        
        C=[C,cc+u >= sum(sum(sum((y1+y).*eo)))+sum(sum(z.*eh))+sum(sum(sum((x1+x).*ev))) + b*sum(sum(sum((x1+x),1),3).*dis)];
        
        optimize(C,obj)
        
        f(n)=obj;
    end
    l=round(lambda/0.2);
    g(l)=mean(f);
    h(l)=std(f);
end
value(h)

[Dyuman Joshi]

What seems to be the problem? If you are getting an error, copy-paste the full error message.

[Qiandong Dong]

There is no error message, and the Matlab solver shows that it is always busy.

[Dyuman Joshi]

I do not have the YALMIP toolbox right now to run your code.

However, you can certainly improve your code -

All the variables defined before x (first use of sdpvar) are constants, bring them out of the loops. (I do not know whether sdpvar variables are constants are not)

And it's not clear to me what is the purpose of the 2nd for loop (n as the loop variable), when there is no change in the code w.r.t n.