📄 fm_opfsdr.m
字号:
mu_Iijcmax= mu(idx+nL); idx = idx + Line.n;
mu_Ijimax = mu(idx+nL); idx = idx + Line.n;
mu_Ijicmax= mu(idx+nL); idx = idx + Line.n;
if Rsrv.n
mu_Prmin = mu(idx+nR); idx = idx + Rsrv.n;
mu_Prmax = mu(idx+nR); idx = idx + Rsrv.n;
mu_sumPrd = mu(idx+1);
end
% Computations for the System: f(thetac,Vc,Qgc,Ps,Pd,lc) = 0
% =====================================================
V_snap = DAE.V; ang_snap = DAE.a;
DAE.V = Vc; DAE.a = angc;
DAE.Gl = zeros(n2,1);
DAE.Gk = zeros(n2,1);
if OPF.line, Line.Y = Y_cont; Line.con = line_cont; end
fm_lf(1);
glambda(SW,1+lc,kg);
glambda(PV,1+lc,kg);
glambda(PQ,1+lc);
Glcall(SW);
Gkcall(SW);
Glcall(PV);
Gkcall(PV);
Glcall(PQ);
% generator reactive powers
DAE.gq = DAE.gq - sparse(busG,1,Qgc,Bus.n,1);
% Demand & Supply
DAE.gp = DAE.gp + sparse(Demand.bus,1,(1+lc)*Pd,Bus.n,1);
DAE.gq = DAE.gq + sparse(Demand.bus,1,(1+lc)*Pd.*qonp,Bus.n,1);
DAE.Gl = DAE.Gl + sparse(Demand.bus,1,Pd,2*Bus.n,1);
DAE.Gl = DAE.Gl + sparse(Demand.bus+Bus.n,1,Pd.*qonp,2*Bus.n,1);
g_Pdc = sparse(Demand.bus,nD,1+lc,2*Bus.n,Demand.n);
g_Pdc = sparse(Demand.bus+Bus.n,nD,(1+lc)*qonp,2*Bus.n,Demand.n);
DAE.gp = DAE.gp - sparse(Supply.bus,1,(1+lc+kg*ksu).*Ps,Bus.n,1);
g_Psc = sparse(Supply.bus,nS,-(1+lc+kg*ksu),2*Bus.n,Supply.n);
DAE.Gl = DAE.Gl - sparse(Supply.bus,1,Ps,2*Bus.n,1);
DAE.Gk = DAE.Gk - sparse(Supply.bus,1,ksu.*Ps,2*Bus.n,1);
gc1p = DAE.gp; gc2p = DAE.glfp;
gc1q = DAE.gq; gc2q = DAE.glfq;
fm_lf(2);
Jlfvc = [DAE.J11, DAE.J12; DAE.J21, DAE.J22];
[Iijc, Jijc, Hijc, Ijic, Jjic, Hjic] = fm_flows(OPF.flow,mu_Iijcmax, mu_Ijicmax);
% Computations for the System: f(theta,V,Qg,Ps,Pd) = 0
% =====================================================
DAE.V = V_snap;
DAE.a = ang_snap;
if OPF.line, Line.Y = Y_orig; Line.con = line_orig; end
fm_lf(1);
gcall(PQ);
glambda(SW,1,0);
glambda(PV,1,0);
% Demand & Supply
DAE.gp = DAE.gp + sparse(Demand.bus,1,Pd,Bus.n,1);
DAE.gq = DAE.gq + sparse(Demand.bus,1,Pd.*qonp,Bus.n,1);
DAE.gp = DAE.gp - sparse(Supply.bus,1,Ps,Bus.n,1);
DAE.gq = DAE.gq - sparse(busG,1,Qg,Bus.n,1);
fm_lf(2);
Jlfv = [DAE.J11, DAE.J12; DAE.J21, DAE.J22];
[Iij, Jij, Hij, Iji, Jji, Hji] = fm_flows(OPF.flow,mu_Iijmax, mu_Ijimax);
% Gradient of [s] variables
% =====================================================
gs = s.*mu - ms;
% Gradient of [mu] variables
% =====================================================
gmu = [Psmin-Ps;Ps-Psmax;Pdmin-Pd;Pd-Pdmax;Qgmin-Qg;Qg-Qgmax; ...
Vmin-DAE.V;DAE.V-Vmax;Qgmin-Qgc;Qgc-Qgmax;Vcmin-Vc;Vc-Vcmax; ...
lcmin-lc;lc-lcmax;Iij-Iijmax;Iijc-Iijcmax;Iji-Iijmax; ...
Ijic-Iijcmax];
if Rsrv.n
gmu(Supply.n+supR) = gmu(Supply.n+supR) + Pr;
gmu = [gmu; Prmin-Pr;Pr-Prmax;sum(Pr)-sum(Pd)];
end
gmu = gmu + s;
% Gradient of [y] = [theta; DAE.V; Qg; Ps; Pd] variables
% =====================================================
if Rsrv.n,
Jg = [Jlfv,g_Qg,g_Ps,g_Pd,Jz1; ...
Jz2,g_Psc,g_Pdc,Jlfvc,DAE.Gk,g_Qgc,DAE.Gl,g_Pr];
else,
Jg = [Jlfv,g_Qg,g_Ps,g_Pd,Jz1; ...
Jz2,g_Psc,g_Pdc,Jlfvc,DAE.Gk,g_Qgc,DAE.Gl];
end
dF_dy = Jg'*ro;
dG_dy(n_d) = -w; % max loading factor
dG_dy(n2+n_gen+1:n2+n_gen+Supply.n) = (1-w)*(Csb + 2*Csc.*Ps + 2*KTBS.*Ps);
dG_dy(n2+n_gen+Supply.n+1:n2+n_a) = -(1-w)*(Cdb + 2*Cdc.*Pd + 2*KTBD.*Pd ...
+ qonp.*(Ddb + 2*qonp.*Ddc.*Pd));
dG_dy(n2+nS) = (1-w)*(Dsb + 2*Dsc.*Qg(nS));
dH_dtV = Jij'*mu_Iijmax + Jji'*mu_Ijimax + [mu_t; mu_Vmax - mu_Vmin];
dH_dtVc = Jijc'*mu_Iijcmax + Jjic'*mu_Ijicmax + [mu_t; mu_Vcmax - mu_Vcmin];
if Rsrv.n,
dH_dy = [dH_dtV; mu_Qgmax-mu_Qgmin; mu_Psmax-mu_Psmin; ...
mu_Pdmax-mu_Pdmin-mu_sumPrd; ...
dH_dtVc; 0; mu_Qgcmax - mu_Qgcmin; mu_lcmax - mu_lcmin; ...
mu_Psmax+mu_Prmax-mu_Prmin+mu_sumPrd];
else
dH_dy = [dH_dtV; mu_Qgmax-mu_Qgmin; mu_Psmax-mu_Psmin; mu_Pdmax-mu_Pdmin; ...
dH_dtVc; 0; mu_Qgcmax - mu_Qgcmin; mu_lcmax - mu_lcmin];
end
gy = dG_dy - dF_dy + dH_dy;
Jh(n_b+1:n_b+Line.n,1:n2) = Jij;
Jh(n_b+1+Line.n:n_b+2*Line.n,1+n2+n_a:n4+n_a) = Jijc;
Jh(n_b+1+2*Line.n:n_b+3*Line.n,1:n2) = Jji;
Jh(n_b+1+3*Line.n:n_b+4*Line.n,1+n2+n_a:n4+n_a) = Jjic;
% Hessian Matrix [D2xLms]
% --------------------------------------------------------------------
H3 = sparse(n_a,n_y);
Hx = -ro(n2+Supply.bus);
Sidx = 1:Supply.n;
H3(n_gen+Sidx,n_d) = Hx;
H3(n_gen+Sidx,n4+1+n_a) = Hx.*ksu;
H5(n_gen+1,n2+n_gen+Sidx) = Hx';
H41(end,n2+n_gen+Sidx) = Hx'.*ksu';
Didx = 1:Demand.n;
Hx = ro(n2+Demand.bus) + qonp.*ro(n3+Demand.bus);
H3(n_gen+Supply.n+Didx,n_d) = Hx;
H5(n_gen+1,n2+n_gen+Supply.n+Didx) = Hx';
H3 = H3 - sparse(n_gen+nS,n2+n_gen+nS,(1-w)*(2*Csc+2*KTBS),n_a,n_y);
H3 = H3 - sparse(nS,n2+nS,(1-w)*2*Dsc,n_a,n_y);
H3 = H3 + sparse(n_gen+Supply.n+nD,n_gen+Supply.n+n2+nD, ...
(1-w)*(2*Cdc+2*KTBD+2*Ddc.*qonp.*qonp),n_a,n_y);
V_snap = DAE.V; ang_snap = DAE.a;
DAE.V = Vc; DAE.a = angc;
if OPF.line, Line.Y = Y_cont; end
Hess_c = -fm_hessian(ro(n2+1:n4))+Hijc+Hjic;
DAE.V = V_snap; DAE.a = ang_snap;
if OPF.line, Line.Y = Y_orig; end
Hess = -fm_hessian(ro(1:n2))+Hij+Hji;
D2xLms = [Hess, H31; -H3; -H41, [Hess_c, Hcolk; Hrowk], H42; -H5];
switch OPF.method
case 1 % Newton Directions
% reduced system
H_m = sparse(1:n_s,1:n_s,mu./s,n_s,n_s);
H_s = sparse(1:n_s,1:n_s,1./s,n_s,n_s);
Jh(:,SW.bus+n2+n_a) = 0;
Jh(:,SW.bus) = 0;
gy = gy+(Jh.')*(H_m*gmu-H_s*gs);
Jd = [D2xLms+(Jh.')*(H_m*Jh),-Jg.';-Jg,Z3];
% reference angle for the actual system
Jd(SW.bus,:) = 0;
Jd(:,SW.bus) = 0;
Jd(SW.bus,SW.bus) = speye(SW.n);
gy(SW.bus) = 0;
% reference angle for the critical system
Jd(:,SW.bus+n2+n_a) = 0;
Jd(SW.bus+n2+n_a,:) = 0;
Jd(SW.bus+n2+n_a,SW.bus+n2+n_a) = speye(SW.n);
gy(SW.bus+n2+n_a) = 0;
% variable increments
Dx = -Jd\[gy; -DAE.gp; -DAE.gq; -gc1p; -gc1q];
Ds = -(gmu+Jh*Dx([1:n_y]));
Dm = -H_s*gs-H_m*Ds;
case 2 % Mehrotra's Predictor-Corrector
% -------------------
% Predictor step
% -------------------
% reduced system
H_m = sparse(1:n_s,1:n_s,mu./s,n_s,n_s);
Jh(:,SW.bus+n2+n_a) = 0;
Jh(:,SW.bus) = 0;
gx = gy+(Jh.')*(H_m*gmu-mu);
Jd = [D2xLms+(Jh.')*(H_m*Jh),-Jg.';-Jg,Z3];
% reference angle for the actual system
Jd(SW.bus,:) = 0;
Jd(:,SW.bus) = 0;
Jd(SW.bus,SW.bus) = speye(SW.n);
gx(SW.bus) = 0;
% reference angle for the critical system
Jd(:,SW.bus+n2+n_a) = 0;
Jd(SW.bus+n2+n_a,:) = 0;
Jd(SW.bus+n2+n_a,SW.bus+n2+n_a) = speye(SW.n);
gx(SW.bus+n2+n_a) = 0;
% LU factorization
[L,U,P] = lu(Jd);
% variable increments
Dx = -U\(L\(P*[gx; -DAE.gp; -DAE.gq; -gc1p; -gc1q]));
Ds = -(gmu+Jh*Dx([1:n_y]));
Dm = -mu-H_m*Ds;
% centering correction
a1 = find(Ds < 0);
a2 = find(Dm < 0);
if isempty(a1), ratio1 = 1; else, ratio1 = -s(a1)./Ds(a1); end
if isempty(a2), ratio2 = 1; else, ratio2 = -mu(a2)./Dm(a2); end
alpha_P = min(1,gamma*min(ratio1));
alpha_D = min(1,gamma*min(ratio2));
c_gap_af = [s + alpha_P*Ds]'*[mu + alpha_D*Dm];
c_gap = s'*mu;
ms = min((c_gap_af/c_gap)^2,0.2)*c_gap_af/n_s;
gs = mu+(Ds.*Dm-ms)./s;
% -------------------
% Corrector Step
% -------------------
% new increment for variable y
gx = gy+(Jh.')*(H_m*gmu-gs);
gx(SW.bus) = 0;
gx(SW.bus+n2+n_a) = 0;
% variable increments
Dx = -U\(L\(P*[gx; -DAE.gp; -DAE.gq; -gc1p; -gc1q]));
Ds = -(gmu+Jh*Dx([1:n_y]));
Dm = -gs-H_m*Ds;
end
% =======================================================================
% Variable Increments
% =======================================================================
Dtheta = Dx(nB); idx = Bus.n; % curr. sys.
DV = Dx(idx+nB); idx = idx + Bus.n;
DQg = Dx(idx+nG); idx = idx + n_gen;
DPs = Dx(idx+nS); idx = idx + Supply.n;
DPd = Dx(idx+nD); idx = idx + Demand.n;
Dthetac = Dx(idx+nB); idx = idx + Bus.n; % crit. sys.
DVc = Dx(idx+nB); idx = idx + Bus.n;
Dkg = Dx(1+idx); idx = idx + 1;
DQgc = Dx(idx+nG); idx = idx + n_gen;
Dlc = Dx(1+idx); idx = idx + 1;
if Rsrv.n, DPr = Dx(idx+nR); idx = idx + Rsrv.n; end
Dro = Dx(1+idx:end); % Lag. mult.
% =======================================================================
% Updating the Variables
% =======================================================================
% Step Lenght Parameters [alpha_P & alpha_D]
%________________________________________________________________________
a1 = find(Ds < 0);
a2 = find(Dm < 0);
if isempty(a1), ratio1 = 1; else, ratio1 = (-s(a1)./Ds(a1)); end
if isempty(a2), ratio2 = 1; else, ratio2 = (-mu(a2)./Dm(a2)); end
alpha_P = min(1,gamma*min(ratio1));
alpha_D = min(1,gamma*min(ratio2));
% New primal variables
%________________________________________________________________________
DAE.a = DAE.a + alpha_P * Dtheta;
DAE.V = DAE.V + alpha_P * DV;
Ps = Ps + alpha_P * DPs;
Pd = Pd + alpha_P * DPd;
Qg = Qg + alpha_P * DQg;
if Rsrv.n, Pr = Pr + alpha_P * DPr; end
angc = angc + alpha_P * Dthetac;
Vc = Vc + alpha_P * DVc;
kg = kg + alpha_P * Dkg;
Qgc = Qgc + alpha_P * DQgc;
lc = lc + alpha_P * Dlc;
s = s + alpha_P * Ds;
% New dual variables
%________________________________________________________________________
ro = ro + alpha_D*Dro;
mu = mu + alpha_D*Dm;
% Objective Function
%________________________________________________________________________
s(find(s == 0)) = epsilon_mu;
Fixd_c = sum(Csa) - sum(Cda) + sum(Dsa) - sum(Dda);
Prop_c = Csb'*Ps - Cdb'*Pd + Dsb'*Qg(nS) - Ddb'*(qonp.*Pd);
TieBreaking = (sum(KTBS.*Ps.*Ps) - sum(KTBD.*Pd.*Pd));
Quad_c = Csc'*(Ps.*Ps) - Cdc'*(Pd.*Pd) - Ddc'*(qonp.*qonp.*Pd.*Pd);
Quad_q = Dsc'*(Qg(nS).*Qg(nS));
if Rsrv.n, Reserve = Cr'*Pr; else, Reserve = 0; end
G_obj = (1-w)*(Fixd_c + Prop_c + Quad_c + Quad_q + TieBreaking + Reserve) - ...
ms*sum(log(s)) - w*lc;
% =======================================================================
% Reducing the Barrier Parameter
% =======================================================================
sigma = max(0.99*sigma, 0.1); % Centering Parameter
c_gap = s'*mu; % Complementarity Gap
ms = min(abs(sigma*c_gap/n_s),1); % Evaluation of the Barrier Parameter
% =======================================================================
% Testing for Convergence
% =======================================================================
test1 = ms <= epsilon_mu;
norma2 = norm(Dx,inf);
test2 = norma2 <= epsilon_2;
norma3 = norm([DAE.gp; DAE.gq; gc1p; gc1q],inf);
test3 = norma3 <= epsilon_1;
norma4 = abs(G_obj-G_obj_k_1)/(1+abs(G_obj));
test4 = norma4 <= epsilon_2;
if test1 & test2 & test3 & test4, break, end
% Displaying Convergence Tests
%________________________________________________________________________
iteration = iteration + 1;
if OPF.show
fm_disp(['Iter. =',fvar(iteration,5),' mu =', fvar(ms,8), ...
' |dy| =', fvar(norma2,8), ' |f(y)| =', ...
fvar(norma3,8),' |dG(y)| =' fvar(norma4,8)])
end
fm_status('opf','update',[iteration, ms, norma2, norma3, norma4], ...
iteration)
if iteration > iter_max, break, end
end
% Some settings ...
%____________________________________________________________________________
warning('on');
Demand = pset(Demand,Pd);
Supply = pset(Supply,Ps);
if Rsrv.n, Rsrv.con(:,10) = Pr; end
Iij = sqrt(Iij);
Iji = sqrt(Iji);
Iijc = sqrt(Iijc);
Ijic = sqrt(Ijic);
Iijmax = sqrt(Iijmax);
Iijcmax = sqrt(Iijcmax);
MVA = Settings.mva;
Pay = ro(1:Bus.n).*DAE.glfp*MVA;
ISOPay = sum(Pay);
% Nodal Congestion Prices (NCPs)
%____________________________________________________________________________
Jlfv(SW.bus,:) = [];
Jlfv(:,SW.bus) = [];
dH_dtV(SW.bus,:) = [];
NCP = -Jlfv'\dH_dtV;
OPF.NCP = [NCP(1:SW.bus-1);0;NCP(SW.bus:end)];
OPF.obj = G_obj;
OPF.ms = ms;
OPF.dy = norma2;
OPF.dF = norma3;
OPF.dG = norma4;
OPF.iter = iteration;
OPF.gpc = gc2p;
OPF.gqc = gc2q;
SNB.init = 0;
LIB.init = 0;
CPF.init = 0;
OPF.init = 2;
% set Pg, Qg, Pl and Ql
Bus.Pl = OPF.basepl*Snapshot(1).Pl + sparse(Demand.bus,1,Pd,Bus.n,1);
Bus.Ql = OPF.basepl*Snapshot(1).Ql + sparse(Demand.bus,1,Pd.*qonp,Bus.n,1);
Bus.Pg = OPF.basepg*Snapshot(1).Pg + sparse(Supply.bus,1,Ps,Bus.n,1);
Bus.Qg = OPF.basepg*Snapshot(1).Qg + sparse(busG,1,Qg,Bus.n,1);
% Display Results
% --------------------------------------------------
TPQ = totp(PQ);
if (Settings.showlf | OPF.show) & clpsat.showopf
OPF.report = cell(1,1);
OPF.report{1,1} = ['Weighting Factor = ',fvar(w,8)];
OPF.report{2,1} = ['Lambda = ',fvar(lc,8)];
OPF.report{3,1} = ['Kg = ',fvar(kg,8)];
OPF.report{4,1} = ['Total Losses = ',fvar(sum(DAE.glfp),8),' [p.u.]'];
OPF.report{5,1} = ['Bid Losses = ',fvar(sum(DAE.glfp)-Snapshot(1).Ploss,8),' [p.u.]'];
if ~noDem
OPF.report{6,1} = ['Total demand = ', ...
fvar(sum(Pd),8),' [p.u.]'];
end
OPF.report{6+(~noDem),1} = ['TTL = ',fvar(sum(Pd)+TPQ,8),' [p.u.]'];
fm_disp
fm_disp('----------------------------------------------------------------')
Settings.lftime = toc;
if Fig.stat, fm_stat; end
if iteration > iter_max
fm_disp('IPM-OPF: Method did not Converge',2)
elseif Fig.main
if ~get(Fig.main,'UserData')
fm_disp('IPM-OPF: Interrupted',2)
else
fm_disp(['IPM-OPF completed in ',num2str(toc),' s'],1)
end
else
fm_disp(['IPM-OPF completed in ',num2str(toc),' s'],1)
if Settings.showlf == 1
fm_stat(OPF.report);
else
if Settings.beep
beep
end
end
end
fm_status('opf','close')
else
if iteration > iter_max
fm_disp('IPM-OPF: Method did not Converge',2)
elseif Fig.main
if ~get(Fig.main,'UserData')
fm_disp(['IPM-OPF: Interrupted'],2)
else
fm_disp(['IPM-OPF completed in ',num2str(toc),' s'],1)
end
else
fm_disp(['IPM-OPF completed in ',num2str(toc),' s'],1)
end
end
if iteration > iter_max,
OPF.conv = 0;
else,
OPF.conv = 1;
end
if Rsrv.n,
OPF.guess = [s; mu; DAE.a; DAE.V; Qg; Ps; Pd; ...
angc; Vc; kg; Qgc; lc; Pr; ro];
else
OPF.guess = [s; mu; DAE.a; DAE.V; Qg; Ps; Pd; ...
angc; Vc; kg; Qgc; lc; ro];
end
OPF.atc = (1+lc)*(sum(Pd)+TPQ)*MVA;
OPF.Vc = Vc;
OPF.ac = angc;
if noDem, Demand = restore(Demand); end
PQ = resetvlim(PQ);
if ~OPF.basepl
PQ = pqreset(PQ,'all');
end
if ~OPF.basepg
SW = swreset(SW,'all');
PV = pvreset(PV,'all');
end
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -