Need help in increasing computational time of quadprog

1 view (last 30 days)
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;
  2 Comments
djibeyrou ba
djibeyrou ba on 15 Dec 2021
I want to run the code. Where i can i download the network data? Please I need your support.
Walter Roberson
Walter Roberson on 15 Dec 2021
If you click on the google URL in the Question, it will allow you to request access.
However, as the user posted the question 10 years ago, they might not want to be bothered to provide the access.

Sign in to comment.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!