Need help in increasing computational time of quadprog

Question

Samir 2011-12-23

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/24637-need-help-in-increasing-computational-time-of-quadprog

评论： Walter Roberson 2021-12-15

I have written a sparse code for optimization using quadprog of the optimization toolbox. It is taking around half an hour for the computational time if number of unknowns are increased. Please people suggest me how can I reduce the simulation time, it could be changing or rewriting anything in the program code like if could have used function here as I am not an expert in Matlab(beginner). Below is the code attached and the excel file used to call data can be downloaded from https://docs.google.com/viewer?a=v&pid=explorer&chrome=true&srcid=0B3w3oa_k44mJODQ1MGE2YzMtNDhjMS00ODVlLWIxMjctM2ViMTA5ZWVkZTA1&hl=en_US

clear;
clc;
tic;
%%Main system data
nb=751;
nl=1011;
h=4;
%Storage and gen assumed at every bus
ng=nb;
nst=nb;
%_______________________________________________
%%Read line data
%Read entire table of line data
B = xlsread('AZ.xlsx','Line Records','b1:b10000');
D = xlsread('AZ.xlsx','Line Records','d1:d10000');
J = xlsread('AZ.xlsx','Line Records','j1:j10000');
%__________________________________________________
%%Formation of Aeq 
aeq=sparse(nb*h+nl*h+1,(3*nb-1+nl)*h);
%__________________________________________________
%Formation of bus dictionary
busdict=zeros(nb,1);
busdict(1)=10435;
running=1;
for iline=1:nl;
    ifrom = B(iline);
    ito = D(iline);
    xline=J(iline);
    ifound=0;
    for i=1:running;
        if busdict(i)==ifrom;
            ifound=i;
            ifromn=i;
        end;
    end;
    if ifound == 0;
        running=running+1;
        busdict(running)=ifrom;
        ifromn=running;
    end;
    ifound=0;
    for i=1:running;
        if busdict(i)==ito;
            ifound=i;
            iton=i;
        end;
    end;
    if ifound == 0
        running = running+1;
        busdict(running)=ito;
        iton=running;
    end;
end;
busdict1=sort(busdict);
%________________________________________________________________________%
inb=sparse(eye(nb*h,nb*h));
aeq(1:nb*h,1:nb*h)=inb;
% %__________________________________________________________________
% %Storage power - assumed at all buses. Zero storage accomodated at LB n UB
aeq(1:nb*h,(nb)*h+1:(nb+nst)*h)=-inb;
% %_____________________________________________________________________
% %Stored Energy at the end of the day should be zero
aeq(nb*h+nl*h+1,(nb)*h+1:(nb+nst)*h)=24/h;
% %______________________________________________________________________
%Line injection power
inh=sparse(eye(h,h));
for k=1:nl;
   stb=B(k);
   stbb=0;
   for look=1:nb;
        if busdict1(look)==stb;
            stbb=look;
        end;
   end;
   endb=D(k);
   endbb=0;
   for look=1:nb;
       if busdict1(look)==endb;
           endbb=look;
       end;
   end;
   % %______________________________________________________________________
% %Line goes from stb to endb
aeq(h*(stbb-1)+1:h*(stbb-1)+h,nb*h+nst*h+1+h*(k-1):nb*h+nst*h+h*(k-1)+h)=-inh;
aeq(h*(endbb-1)+1:h*(endbb-1)+h,nb*h+nst*h+1+h*(k-1):nb*h+nst*h+h*(k-1)+h)=inh;
end;
for k=1:nl;
    x=J(k);
    stb=B(k);
    stbb=0;
    for look=1:nb;
        if busdict1(look)==stb;
            stbb=look;
        end;
    end;
    endb=D(k);
    endbb=0;
    for look=1:nb;
        if busdict1(look)==endb;
            endbb=look;
        end
    end
    if stbb~=1;
       aeq(nb*h+1+(k-1)*h:nb*h+(k-1)*h+h,(2*nb*h+nl*h)+1+(stbb-2)*h:(2*nb*h+nl*h)+(stbb-2)*h+h)=-1/x*inh;
    end
    if endbb~=1;
       aeq(nb*h+1+(k-1)*h:nb*h+(k-1)*h+h,(2*nb*h+nl*h)+1+(endbb-2)*h:(2*nb*h+nl*h)+(endbb-2)*h+h)=1/x*inh;
    end
end
% % %______________________________________________________________________
% %Power Flow vs. delta - this completes the formation of matrix aeq
inl=sparse(eye(nl*h,nl*h));
aeq(nb*h+1:(nb+nl)*h,nb*h+nst*h+1:(nb+nst+nl)*h)=inl;
% % %_______________________________________________________________________
% %%Formation of beq
beq=zeros(nb*h+nl*h+1,1);
U=xlsread('AZ.xlsx','Bus Records','n2:n10000');%Load 
V=xlsread('AZ.xlsx','Bus Records','o2:o10000');%Wind
for k=1:nb*h
    load=U(k)/100;
    wind=V(k)/100;
    beq(k,1)=load-wind;
end
%%Formation of A
a=sparse(nl*h*2+nst*h,ng*h+nl*h+nst*h+(nb-1)*h);
a(1:nl*h,(ng+nst)*h+1:(ng+nst+nl)*h)=inl;
a(nl*h+1:(nl+nl)*h,(ng+nst)*h+1:(ng+nst+nl)*h)=-inl;
n=0;m=0;
for i=1:nst;
a(2*nl*h+1+m:2*nl*h+(h-1)+m,ng*h+1+n:ng*h+h-1+n)=sparse(tril(ones(h-1,h-1)))*24/h;
k=a(2*nl*h+1+m:2*nl*h+(h-1)+m,ng*h+1+n:ng*h+h-1+n);
    if h==3;
        a(2*nl*h+h+m:2*nl*h+h+m,ng*h+1+n:ng*h+h-1+n)=-k(2:h-1,1:h-1);
        j=m;
    else 
        a(2*nl*h+h+m:2*nl*h+h+h-3+m,ng*h+1+n:ng*h+h-1+n)=-k(2:h-1,1:h-1);
        j=m;
    end
n=n+h;m=2*h+m-3;
end
%Maximum stored power
inb=sparse(eye(nb*h,nb*h));
a(2*nl*h+2*h-2+j:2*nl*h+2*h-3+j+nst*h,(nb)*h+1:(nb+nst)*h)=inb;
b=zeros(nl*h*2+nst*h,1);
%Reading Line constarints from excel
Y=xlsread('AZ.xlsx','Line Records','l1:l10000');
k=0;
for j=1:nl;
    kline=Y(j)/100;
    b(1+k:h+k,1)=kline;
    k=k+h;
end
k=0;
for j=1:nl;
    kline=Y(j)/100;
    b(nl*h+1+k:nl*h+h+k,1)=kline;
    k=k+h;
end
%Energy stored
S=xlsread('AZ.xlsx','Storage','q1:q10000');
m=0;j=0;n=0;
for i=1:nst;
    s=S(i)/100;
    b(2*nl*h+1+m:2*nl*h+h-1+m,1)=ones(h-1,1)*s;
    if h==3;
        b(2*nl*h+h+m:2*nl*h+h+m,1)=ones(h-2,1)*0;
        j=m;
        %n=n+h;m=2*h+m-3;
    else 
        b(2*nl*h+h+m:2*nl*h+h+h-3+m,1)=ones(h-2,1)*0;
        j=m;
        %n=n+h;m=2*h+m-3;
    end
m=2*h+m-3;
end
%Maximum Power stored
P=xlsread('AZ.xlsx','Storage','p1:p10000');
k=0;
for i=1:nst;
    p=P(i)/100;
    b(2*nl*h+2*h-2+j+k:2*nl*h+3*h-3+j+k,1)=ones(h,1)*p;
    k=k+h;
end
%%Cost function to be minimized
%c=zeros(1,ng*h+nl*h+nst*h+(nb-1)*h);
H=zeros(ng*h+nl*h+nst*h+(nb-1)*h,ng*h+nl*h+nst*h+(nb-1)*h);
Q=xlsread('AZ.xlsx','Storage','s2:s10000');
C=xlsread('AZ.xlsx','Bus Records','s2:s10000');
n=0;k=0;
for ro=1:h:ng*h;
    q=Q(ro+k-n);
    for added=1:h;
        H(ro+added-1,ro+added-1)=2*q*100;
    end
    k=k+1;
    n=n+h;
end
f=C*100;
for i=ng*h+1:ng*h+nl*h+nst*h+(nb-1)*h;
    f(i,1)=0;
end
% for g=1:nb*h;
%     lmp=C(g);
%     c(1,g)=lmp*(24/h);
% end;
%%Construct lb and ub vectors
lb=zeros((3*nb-1+nl)*h,1);
esratg=xlsread('AZ.xlsx','Storage','p1:p10000');
k=0;n=0;
for ro=nb*h+1:h:2*nb*h;
    lb(ro)=0;
    es=ro-nb*h-h*k+n;
    %es=((ro-h*(n+2)-1)/h+k);
    rate=esratg(es);
    for added=1:h-1;
        lb(ro+added)=-rate;
    end;
    k=k+1;
    n=n+1;
end;
n=0;k=0;
for ro = 2*nb*h+1:h:2*nb*h+nl*h;
    y=-Y(ro-2*nb*h-h*k+n);%-Y((ro-h*(n+5)-1)/h+k);
    for added=1:h;
    lb(ro+added-1)=y;
    end
    k=k+1;
    n=n+1;
end
lb(2*nb*h+nl*h+1:(3*nb-1+nl)*h,1)=-pi/6;    
ub=-lb;
G=xlsread('AZ.xlsx','Storage','r2:r10000');
n=0;k=0;
for ro=1:h:nb*h;
    g=G(ro+k-n);
    for added=1:h;
        ub(ro+added-1)=g;
    end
    k=k+1;
    n=n+h;
end
k=0;
n=0;
for ro=nb*h+1:h:2*nb*h;
    ub(ro+h-1)=0;
    es=ro-nb*h-h*k+n;
    %es=((ro-h*(n+2)-1)/h+k);
    rate=esratg(es);
    for added=1:h-1;
        ub(ro+added-1)=rate;
    end;
    k=k+1;
    n=n+1;
end;
%ub(1:982)=Inf;
%aeq(1479:1488,:)=0;beq(1479:1488)=0;
% aeq(3433:3434,:)=0;
% beq(3433:3434)=0;
% aeq(1499+139*2:1498+140*2,:)=0;beq(1499+139*2:1498+140*2)=0;
% aeq(1499+182*2:1498+183*2,:)=0;beq(1499+182*2:1498+183*2)=0;
options=optimset('Algorithm','interior-point-convex');
[X,fval]=quadprog(H,f,a,b,aeq,beq,lb,ub,[],options);
%Aeq=full(aeq);
% toc;