Need help in increasing computational time of quadprog
조회 수: 3 (최근 30일)
이전 댓글 표시
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
2021년 12월 15일
I want to run the code. Where i can i download the network data? Please I need your support.
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 Center 및 File 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!