📄 fm_gams.m
字号:
[Piji,Pjii,Qgi] = updatePF(Pd(i,:)',Ps(i,:)',iBQg);
a(i,:) = DAE.a';
V(i,:) = DAE.V';
Pij(i,:) = Piji';
Pji(i,:) = Pjii';
Qg(i,iBQg) = Qgi';
end
ro = MCP;
end
% ------------------------------------------------------------------
case '3' % S T A N D A R D O P T I M A L P O W E R F L O W
% ------------------------------------------------------------------
if type == 1 % Single Period Auction
[Ps,Pd,V,a,Qg,ro,Pij,Pji,mV,mFij,mFji,modelstat,solvestat] = ...
psatgams(file,nBus,nLine,nPs,nPd,nSW,nPV,nBusref,control, ...
flow,Gh,Bh,Li,Lj,Ps_idx,Pd_idx,S,D,X,L);
NCP = compNCP(V,a,mV,mFij,mFji);
elseif ~rem(type,2) % Single/Multi Period Auction with UC
[Ps,Pd,V,a,Qg,ro,Pij,Pji,mV,mFij,mFji,modelstat,solvestat] = ...
psatgams(file,nBus,nLine,nPs,nPd,nSW,nPV,nBusref,nH,control, ...
flow,Gh,Bh,Li,Lj,Ps_idx,Pd_idx,S,D,X,L,Ch);
NCP = zeros(length(Ps(:,1)),Bus.n);
for i = 1:length(Ps(:,1))
NCPi = compNCP(V(i,:)',a(i,:)',mV(i,:)',mFij(i,:)',mFji(i,:)');
NCP(i,:) = NCPi';
end
end
% ------------------------------------------------------------------
case '7' % C O N G E S T I O N M A N A G E M E N T
% ------------------------------------------------------------------
%lambda_values = [0.0:0.01:0.61];
%n_lambda = length(lambda_values);
%GAMS.dpgup = zeros(Supply.n,n_lambda);
%GAMS.dpgdw = zeros(Supply.n,n_lambda);
%GAMS.dpdup = zeros(Demand.n,n_lambda);
%GAMS.dpddw = zeros(Demand.n,n_lambda);
%for i = 1:length(lambda_values)
%lambda.val(1) = lambda_values(i);
iteration = 0;
idx_gen = zeros(Supply.n,1);
Psc_idx = Ps_idx;
Psm_idx = zeros(Bus.n,Supply.n);
while 1
[Ps,Pd,dPSup,dPSdw,dPDup,dPDdw,V,a,Qg,ro,Pij,Pji,mV,mFij,mFji, ...
lambdac,kg,Vc,ac,Qgc,Pijc,Pjic,Pceq,lambdam,modelstat,solvestat] = ...
psatgams('fm_cong',nBus,nLine,nPs,nPd,nBusref,nSW,nPV,control,flow, ...
Gh,Bh,Ghc,Bhc,Li,Lj,Ps_idx,Psc_idx,Psm_idx,Pd_idx, ...
SW_idx,PV_idx,S,D,X,L,lambda);
iteration = iteration + 1;
if iteration > 10
fm_disp('* * * Maximum number of iteration with no convergence!')
break
end
idx = psupper(Supply,(1+lambdac+kg)*Ps);
if sum(idx_gen(idx)) == length(idx)
fm_disp(['* * * iter = ',num2str(iteration), ...
', #viol., ',num2str(length(find(idx_gen))), ...
' lambda = ', num2str(lambdac), ...
' kg = ', num2str(kg)])
break
else
% loop until there are no violations of power supply limits
idx_gen(idx) = 1;
fm_disp(['* * * iter = ',num2str(iteration),', #viol. = ', ...
num2str(length(idx)),', lambda = ', ...
num2str(lambdac),' kg = ', num2str(kg)])
drawnow;
Psc_idx = psidx(Supply,~idx_gen);
Psm_idx = psidx(Supply,idx_gen);
end
end
%GAMS.dpgup(:,i) = dPSup;
%GAMS.dpgdw(:,i) = dPSdw;
%GAMS.dpdup(:,i) = dPDup;
%GAMS.dpddw(:,i) = dPDdw;
%if ~rem(i,10), disp(['Current lambda = ',num2str(lambda_values(i))]),end
%end
%GAMS.lvals = lambda_values;
NCP = compNCP(V,a,mV,mFij,mFji);
% ------------------------------------------------------------------
case '4' % V O L T A G E S T A B I L I T Y C O N S T R A I N E D
% O P T I M A L P O W E R F L O W
% ------------------------------------------------------------------
if type == 1 % Single Period Auction
[Ps,Pd,V,a,Qg,ro,Pij,Pji,mV,mFij,mFji, ...
lambdac,kg,Vc,ac,Qgc,Pijc,Pjic,Pceq,modelstat,solvestat] = ...
psatgams(file,nBus,nLine,nPs,nPd,nBusref,nSW,nPV,control,flow, ...
Gh,Bh,Ghc,Bhc,Li,Lj,Ps_idx,Pd_idx,SW_idx,PV_idx,S,D,X,L,lambda);
NCP = compNCP(V,a,mV,mFij,mFji);
elseif ~rem(type,2) % Single/Multi Period Auction with UC
[Ps,Pd,V,a,Qg,ro,Pij,Pji,mV,mFij,mFji, ...
lambdac,kg,Vc,ac,Qgc,Pijc,Pjic,Pceq,modelstat,solvestat] = ...
psatgams(file,nBus,nLine,nPs,nPd,nBusref,nH,control,flow, ...
Gh,Bh,Ghc,Bhc,Li,Lj,Ps_idx,Pd_idx,S,D,X,L,lambda,Ch);
NCP = zeros(length(lambdac),Bus.n);
for i = 1:length(lambdac)
NCPi = compNCP(V(i,:)',a(i,:)',mV(i,:)',mFij(i,:)',mFji(i,:)');
NCP(i,:) = NCPi';
end
elseif type == 3 % Pareto Set Single Period Auction
fm_disp
for i = 1:length(omega)
fm_disp(sprintf(' VSC-OPF #%d, %3.1f%% - omega: %5.4f', ...
i,100*i/length(omega),omega(i)))
lambda.val = [GAMS.lmin(1),GAMS.lmax(1),omega(i),GAMS.line];
[Psi,Pdi,Vi,ai,Qgi,roi,Piji,Pjii,mV,mFij,mFji, ...
lambdaci,kgi,Vci,aci,Qgci,Pijci,Pjici,Pceq,modelstat,solvestat] = ...
psatgams(file,nBus,nLine,nPs,nPd,nBusref,nSW,nPV,control,flow, ...
Gh,Bh,Ghc,Bhc,Li,Lj,Ps_idx,Pd_idx,SW_idx,PV_idx, ...
S,D,X,L,lambda);
gams_mstat(modelstat)
gams_sstat(solvestat)
Ps(i,:) = Psi';
Pd(i,:) = Pdi';
V(i,:) = Vi';
a(i,:) = ai';
Qg(i,:) = Qgi';
ro(i,:) = roi';
Pij(i,:) = Piji';
Pji(i,:) = Pjii';
lambdac(i,:) = lambdaci';
kg(i,:) = kgi';
Vc(i,:) = Vci';
ac(i,:) = aci';
Qgc(i,:) = Qgci';
Pijc(i,:) = Pijci';
Pjic(i,:) = Pjici';
NCPi = compNCP(Vi,ai,mV,mFij,mFji);
NCP(i,:) = NCPi';
end
fm_disp
end
% ------------------------------------------------------------------
case '5' % M A X I M U M L O A D I N G C O N D I T I O N
% ------------------------------------------------------------------
if type == 1 % Single Period Auction
[Ps,Pd,V,a,Qg,ro,Pij,Pji,lambdac,kg,modelstat,solvestat] = ...
psatgams(file,nBus,nLine,nPs,nPd,nSW,nPV,nBusref, ...
control,flow,Gh,Bh,Li,Lj,Ps_idx,Pd_idx, ...
SW_idx,PV_idx,S,D,X,L);
elseif ~rem(type,2) % Single/Multi Period Auction with UC
[Ps,Pd,V,a,Qg,ro,Pij,Pji,lambdac,kg,modelstat,solvestat] = ...
psatgams(file,nBus,nLine,nPs,nPd,nSW,nPV,nBusref,nH, ...
control,flow,Gh,Bh,Li,Lj,Ps_idx,Pd_idx,SW_idx, ...
PV_idx,S,D,X,L,Ch);
end
% ------------------------------------------------------------------
case '6' % C O N T I N U A T I O N
% O P T I M A L P O W E R F L O W
% ------------------------------------------------------------------
initial_time = clock;
if type == 1 % single period OPF, no discrete variables
% number of steps
i = 0;
last_point = 0;
% initial lambda = 0. Base case has to be feasible
lmin = 0;
lmax = 0;
Lambda = lmin;
% save actual CPF settings
CPF_old = CPF;
%CPF.nump = 50;
CPF.show = 0;
CPF.type = 3;
CPF.sbus = 0;
CPF.vlim = 1;
CPF.ilim = 1;
CPF.qlim = 1;
CPF.init = 0;
CPF.step = 0.25;
control = '6';
% save actual power flow data
% ------------------------------------------------------
snappg = Snapshot(1).Pg;
%Snapshot(1).Pg = [];
J11old = DAE.J11;
J12old = DAE.J12;
J21old = DAE.J21;
J22old = DAE.J22;
Bus_old = Bus;
% defining voltage limits
[Vmax,Vmin] = fm_vlim(1.2,0.8);
fm_disp
stop_opf = 0;
while 1
% OPF step
% ------------------------------------------------------
i = i + 1;
fm_disp(sprintf('Continuation OPF #%d, lambda_c = %5.4f', ...
i,lmin))
lambda.val = [lmin,lmax,0,GAMS.line];
% call GAMS
[Psi,Pdi,Vi,ai,Qgi,roi,Piji,Pjii,mV,mFij,mFji, ...
lambdaci,kgi,Vci,aci,Qgci,Pijci,Pjici,mPceq,ml,modelstat,solvestat] = ...
psatgams(file,nBus,nLine,nPs,nPd,nBusref,nSW,nPV,control,flow, ...
Gh,Bh,Ghc,Bhc,Li,Lj,Ps_idx,Pd_idx,SW_idx,PV_idx,S,D,X,L,lambda);
gams_mstat(modelstat)
gams_sstat(solvestat)
Lambda(i,1) = lambdaci;
Ps(i,:) = Psi';
Pd(i,:) = Pdi';
V(i,:) = Vi';
a(i,:) = ai';
Qg(i,:) = Qgi';
ro(i,:) = roi';
Pij(i,:) = Piji';
Pji(i,:) = Pjii';
lambdac(i,:) = lambdaci;
kg(i,:) = kgi;
Vc(i,:) = Vci';
ac(i,:) = aci';
Qgc(i,:) = Qgci';
Pijc(i,:) = Pijci';
Pjic(i,:) = Pjici';
NCPi = compNCP(Vi,ai,mV,mFij,mFji);
NCP(i,:) = NCPi';
ML(i,1) = ml;
% check consistency of the solution (LMP > 0)
if modelstat > 3 %min(abs(roi)) < 1e-5
fm_disp('Unfeasible OPF solution. Discarding last solution.')
Lambda(end) = [];
Ps(end,:) = [];
Pd(end,:) = [];
V(end,:) = [];
a(end,:) = [];
Qg(end,:) = [];
ro(end,:) = [];
Pij(end,:) = [];
Pji(end,:) = [];
lambdac(end) = [];
kg(end) = [];
Vc(end,:) = [];
ac(end,:) = [];
Qgc(end,:) = [];
Pijc(end,:) = [];
Pjic(end,:) = [];
NCP(end,:) = [];
ML(end) = [];
lambdaci = lambdac(end,:);
break
end
% ------------------------------------------------------
% Bid variations to allow loading parameter increase
%
% d mu_Pceq_i
% D P_i = -sign(P_i) -------------
% d mu_lambda
%
% where:
%
% P_i = power bid i
% mu_Pceq_i = Lagrangian multiplier of critical PF eq. i
% mu_lambda = Lagrangian multiplier of lambda
% ------------------------------------------------------
delta = 0.05;
while 1
if abs(ml) > 1e-5
deltaPd = ml./mPceq(Demand.bus)/(1+lambdaci);
deltaPs = -ml./mPceq(Supply.bus)/(1+lambdaci+kgi);
delta_max = norm([deltaPs; deltaPd]);
if delta_max == 0, delta_max = 1; end
deltaPd = deltaPd/delta_max;
deltaPs = deltaPs/delta_max;
else
deltaPd = zeros(Demand.n,1);
deltaPs = zeros(Supply.n,1);
end
%ml
%mPceq
%delta_max = max(norm([deltaPs; deltaPd]));
%if delta_max == 0, delta_max = 1; end
DPs(i,:) = deltaPs'/Settings.mva;
DPd(i,:) = deltaPd'/Settings.mva;
Pdi = pdbound(Demand,Pd(i,:)' + delta*deltaPd.*Pd(i,:)');
Psi = psbound(Supply,Ps(i,:)' + delta*deltaPs.*Ps(i,:)');
% CPF step
% ------------------------------------------------------
if GAMS.basepl
PQ = pqreset(PQ,'all');
PV = pvreset(PV,'all');
SW = swreset(SW,'all');
Snapshot(1).Pg = snappg;
else
PQ = pqzero(PQ,'all');
PV = pvzero(PV,'all');
SW = swzero(SW,'all');
Snapshot(1).Pg = zeros(Bus.n,1);
end
Demand = pset(Demand,Pdi);
pqsum(Demand,1);
Supply = pset(Supply,Psi);
pgsum(Supply,1);
swsum(Supply,1);
DAE.V = Vci;
DAE.a = aci;
PV = setvg(PV,'all',DAE.V(PV.bus));
SW = setvg(SW,'all',DAE.V(SW.bus));
DAE.x = Snapshot(1).x;
Bus.Pg = Bus_old.Pg;
Bus.Qg = Bus_old.Qg;
Bus.Pl = Bus_old.Pl;
Bus.Ql = Bus_old.Ql;
% avoid aborting CPF routine due to limits
% ------------------------------------------------------
% voltage limits
idx = find(abs(DAE.V-Vmax) < CPF.tolv | DAE.V > Vmax);
if ~isempty(idx)
DAE.V(idx) = Vmax(idx)-1e-6-CPF.tolv;
end
idx = find(abs(DAE.V-Vmin) < CPF.tolv | DAE.V < Vmin);
if ~isempty(idx)
DAE.V(idx) = Vmin(idx)+1e-6+CPF.tolv;
end
CPF.kg = 0;
CPF.lambda = 1; %lambdaci + 1;
CPF.linit = 1+lambdaci*0.25;
CPF.init = 0;
% set contingency for CPF analysis
if GAMS.line
Line.con = line_cont;
fm_y;
end
% ---------------------------------------------
% call continuation power flow routine
fm_cpf('gams');
%CPF.lambda = CPF.lambda + 1;
% ---------------------------------------------
% reset admittance line
if GAMS.line
Line.Y = Y_orig;
Line.con = line_orig;
end
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -