Need help in increasing computational time of quadprog

조회 수: 3 (최근 30일)
Samir
Samir 2011년 12월 23일
댓글: 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;
  댓글 수: 2
djibeyrou ba
djibeyrou ba 2021년 12월 15일
I want to run the code. Where i can i download the network data? Please I need your support.
Walter Roberson
Walter Roberson 2021년 12월 15일
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.

댓글을 달려면 로그인하십시오.

답변 (0개)

카테고리

Help CenterFile Exchange에서 Quadratic Programming and Cone Programming에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by