📄 fm_cpf.m
字号:
DAE.V = DAE.V + inc_corr(DAE.n+Bus.n+1: DAE.n+2*Bus.n);
kg = kg + inc_corr(end-1);
lambda = lambda + inc_corr(end);
iter_corr = iter_corr + 1;
err_corr = max(abs(inc_corr));
end
% Generator reactive power computations
if qlim
DAE.gq = zeros(Bus.n,1);
fm_call('pq');
glambda(Demand,lambda)
Bus.Ql = DAE.gq;
Q = DAE.glfq;
Bus.Qg = Q + Bus.Ql;
[Qmax_idx,Qmin_idx] = pvlim(PV,Bus.Qg(PV.bus));
[Qswmax_idx,Qswmin_idx] = swlim(SW,Bus.Qg(SW.bus));
if ~kqmin
Qmin_idx = [];
Qswmin_idx = [];
end
end
[PQ,lim_v] = pqlim(PQ,vlim,sp,lambda,one);
if lim_v
sigma_corr = 1;
proceed = 1;
if stop > 1
sigma_corr = -1;
break
end
end
if ilim
VV = DAE.V.*exp(jay*DAE.a);
switch ilim
case 1
Fij = abs(((VV(Line.from) - tps.*VV(Line.to)).*y + ...
VV(Line.from).*(jay*chrg))./(tps.*conj(tps)));
Fji = abs((VV(Line.to) - VV(Line.from)./tps).*y + ...
VV(Line.to).*(jay*chrg));
case 2
Fij = abs(real(VV(Line.from).*conj(((VV(Line.from) - ...
tps.*VV(Line.to)).*y + ...
VV(Line.from).*(jay*chrg))./(tps.*conj(tps)))));
Fji = abs(real(VV(Line.to).*conj((VV(Line.to) - ...
VV(Line.from)./tps).*y + VV(Line.to).*(jay*chrg))));
case 3
Fij = abs(VV(Line.from).*conj(((VV(Line.from) - ...
tps.*VV(Line.to)).*y + VV(Line.from).* ...
(jay*chrg))./(tps.*conj(tps))));
Fji = abs(VV(Line.to).*conj((VV(Line.to) - ...
VV(Line.from)./tps).*y + VV(Line.to).*(jay*chrg)));
end
Iij_idx = find(Fij > Imax & Iij);
Iji_idx = find(Fji > Imax & Iji);
end
% determination of the initial loading factor in case
% of infeasible underexcitation of generator at zero
% load condition
if ~iterazione & qlim & ~isempty(Qmin_idx) & count_qmin <= 5
count_qmin = count_qmin+1;
if count_qmin > 5
fm_disp([sp,'There are generator Qmin limit violations at ' ...
'the initial point.'])
fm_disp([sp,'Generator Qmin limits will be disabled.'])
kqmin = 0;
lambda_old = CPF.linit;
lambda = CPF.linit;
sigma_corr = 1;
else
lambda_old = lambda_old + d_lambda;
fm_disp([sp,'Initial loading parameter changed to ', ...
fvar(lambda_old-one,4)])
end
proceed = 0;
break
end
% Check for Hopf Bifurcations
if DAE.n >= 2 & CPF.hopf
As = DAE.Fx-DAE.Fy*(DAE.Jlfv\DAE.Gx);
opt.disp = 0;
auto = eigs(As,2,'SR',opt);
auto = round(auto/Settings.lftol)*Settings.lftol;
hopf_idx = find(real(auto) > 0);
if ~isempty(hopf_idx)
if imag(auto(hopf_idx)) ~= 0
fm_disp([sp,'Hopf bifurcation encountered.'])
end
end
end
if ~isempty(Iij_idx) & ilim
Iij_idx = Iij_idx(1);
Iij(Iij_idx) = 0;
fm_disp([sp,fw,'from bus #',fvar(Line.from(Iij_idx),4), ...
' to bus #',fvar(Line.to(Iij_idx),4), ...
' reached I_max at lambda = ',fvar(lambda-one,9)],1)
sigma_corr = 1;
proceed = 1;
if stop > 1
proceed = 1;
sigma_corr = -1;
break
end
end
if ~isempty(Iji_idx) & ilim
Iji_idx = Iji_idx(1);
Iji(Iji_idx) = 0;
fm_disp([sp,fw,'from bus #',fvar(Line.to(Iji_idx),4), ...
' to bus #',fvar(Line.from(Iji_idx),4), ...
' reached I_max at lambda = ',fvar(lambda-one,9)],1)
sigma_corr = 1;
proceed = 1;
if stop > 1
proceed = 1;
sigma_corr = -1;
break
end
end
if lambda < CPF.linit
cpfmsg = [sp,'lambda is lower than initial value'];
if iterazione > 5
proceed = 0;
break
end
end
if abs(lambda-lambda_old) > 5*abs(d_lambda) & perp & iterazione
fm_disp([sp,'corrector solution is too far from predictor value'])
proceed = 0;
break
end
if lambda > lambda_old & lambda < max(l_vect)
fm_disp([sp,'corrector goes back increasing lambda'])
proceed = 0;
break
end
lim_q = (~isempty(Qmax_idx) | ~isempty(Qmin_idx) | ...
(~isempty(Qswmax_idx) | ~isempty(Qswmin_idx))*type) & qlim;
lim_i = (~isempty(Iij_idx) | ~isempty(Iji_idx)) & ilim;
anylimit = (lim_q | lim_v | lim_i) & CPF.stepcut;
if iter_corr > iter_max % no convergence
if lambda_old < 0.5*max(l_vect)
fm_disp([sp,'Convergence problems in the unstable curve'])
lambda = -1;
proceed = 1;
break
end
if sigma_corr < 1e-3
proceed = 1;
break
end
if CPF.show
fm_disp([sp,'Max # of iters. at corrector step.'])
fm_disp([sp,'Reduced Variable Increments in ', ...
'Corrector Step ', num2str(0.5*sigma_corr)])
end
cpfmsg = [sp,'Reached maximum number of iterations ', ...
'for corrector step.'];
proceed = 0;
break
elseif anylimit & sigma_corr > 0.11
proceed = 0;
break
elseif lim_q
if ~isempty(Qmax_idx)
PQ = add(PQ,pqdata(PV,Qmax_idx(1),'qmax',sp,lambda,one,noSup));
if noSup,
PV = move2sup(PV,Qmax_idx(1));
else
PV = remove(PV,Qmax_idx(1));
end
elseif ~isempty(Qmin_idx)
PQ = add(PQ,pqdata(PV,Qmin_idx(1),'qmin',sp,lambda,one,noSup));
if noSup
PV = move2sup(PV,Qmin_idx(1));
else
PV = remove(PV,Qmin_idx(1));
end
elseif ~isempty(Qswmax_idx)
PQ = add(PQ,pqdata(SW,Qswmax_idx(1),'qmax',sp,lambda,one));
SW = remove(SW,Qswmax_idx(1));
elseif ~isempty(Qswmin_idx)
PQ = add(PQ,pqdata(SW,Qswmin_idx(1),'qmin',sp,lambda,one));
SW = remove(SW,Qswmin_idx(1));
end
Qmax_idx = [];
Qmin_idx = [];
Qswmax_idx = [];
Qswmin_idx = [];
if ~iterazione
d_lambda = 0;
lambda_old = CPF.linit;
lambda = CPF.linit;
if perp
inc(end) = 1e-5;
else
Ljac(end) = 1;
end
else
lambda = lambda_old;
sigma_corr = 1;
proceed = 0;
break
end
else
proceed = 1;
sigma_corr = 1;
if stop & ksign < 0,
fm_disp([sp,'Saddle-node Bifurcation or Limit ', ...
'Induced Bifurcation encountered.'])
sigma_corr = -1;
end
break
end
end
if proceed == 1
if lambda < 0 & iterazione > 1
fm_disp([sp,'lambda < 0 at iteration ',num2str(iterazione)])
break
end
l_vect = [l_vect; lambda];
if CPF.show
fm_disp(['Iteration = ',fvar(iterazione+1,5),'lambda =', ...
fvar(lambda-one,9), ' kg =',fvar(kg,9)],1)
end
iterazione = iterazione + 1;
%etaR(iterazione) = min(diag(DAE.J22-DAE.J21*(DAE.J11\DAE.J12)));
if Syn.n
Syn.pm(syn_idx) = (lambda+kg)*Pm_syn(syn_idx);
end
fm_out(2,lambda,iterazione)
if Syn.n
Syn.pm(syn_idx) = Pm_syn(syn_idx);
end
fm_status('cpf','update',[lambda, DAE.V(PQidx)],iterazione)
if sigma_corr < tol, break, end
sigma_corr = 1;
if lambda > lambda_old
V_crit = DAE.V;
a_crit = DAE.a;
x_crit = DAE.x;
k_crit = kg;
l_crit = lambda;
end
elseif proceed == 0
if abs(lambda-CPF.linit) < 0.001 & ...
abs(lambda-lambda_old) <= abs(d_lambda) & ...
iterazione > 1
fm_disp([sp,'d_lambda = ',num2str([d_lambda, lambda])])
break
end
DAE.V = V_old;
DAE.a = ang_old;
DAE.x = x_old;
kg = kg_old;
lambda = lambda_old;
sigma_corr = 0.1*sigma_corr;
if sigma_corr < tol,
if ~isempty(cpfmsg)
fm_disp(cpfmsg)
end
if iterazione == 0
fm_disp([sp,'Infeasible initial loading condition.'])
else
fm_disp([sp,'Convergence problem encountered.'])
end
break
end
end
% stop routine
% --------------------------------------------------------------
if iterazione >= CPF.nump & ~strcmp(fun,'gams')
fm_disp([sp,'Reached maximum number of points.'])
break
end
if Fig.main
if ~get(Fig.main,'UserData'), break, end
end
% predictor step
% --------------------------------------------------------------
DAE.J11 = sparse(Bus.n,Bus.n);
DAE.J21 = sparse(Bus.n,Bus.n);
DAE.J12 = sparse(Bus.n,Bus.n);
DAE.J22 = sparse(Bus.n,Bus.n);
fm_lf(2);
if noSup
Gyreactive(PV)
else
Gycall(PV);
end
Gyreactive(SW)
DAE.Jlfv = [DAE.J11, DAE.J12; DAE.J21, DAE.J22];
if Bus.n > 250
[L,U,P] = luinc([DAE.Fx,DAE.Fy,Fk;DAE.Gx,DAE.Jlfv,DAE.Gk;kjac],1e-6);
else
[L,U,P] = lu([DAE.Fx,DAE.Fy,Fk;DAE.Gx,DAE.Jlfv,DAE.Gk;kjac]);
end
dz_dl = -U\(L\(P*[Flambda;DAE.Gl;0]));
Jsign_old = Jsign;
Jsign = signperm(P)*sign(prod(diag(U)));
if Jsign ~= Jsign_old & Sflag & iterazione > 1
ksign = -1;
Sflag = 0;
end
Norm = norm(dz_dl,inf);
if Norm == 0, Norm = 1; end
d_lambda = ksign*CPF.step/Norm;
d_lambda = min(d_lambda,0.35);
d_lambda = max(d_lambda,-0.35);
if ksign > 0, d_lambda = max(d_lambda,0); end
if ksign < 0, d_lambda = min(d_lambda,0); end
d_z = d_lambda*dz_dl;
inc = [d_z; d_lambda];
if ~perp & iterazione
a = inc(DAE.n+Bus.n+PQidx);
a = min(a,0.025);
a = max(a,-0.025);
inc(DAE.n+Bus.n+PQidx) = a;
Ljac(end) = 0;
Ljac(DAE.n+Bus.n+PQidx) = 1;
end
end
fm_out(3,0,iterazione)
[lambda_crit, idx_crit] = max(Varout.t);
if isnan(lambda_crit), lambda_crit = lambda; end
if isempty(lambda_crit), lambda_crit = lambda; end
if CPF.show
fm_disp([sp,'Maximum Loading Parameter lambda_max = ', ...
fvar(lambda_crit-one,9)],1)
end
% Visualization of results
% --------------------------------------------------------------
if nodyn
DAE.n = 0;
Varout.idx = Varout.idx-1;
end
Settings.lftime = toc;
Settings.iter = iterazione;
DAE.V = V_crit;
DAE.a = a_crit;
DAE.x = x_crit;
kg = k_crit;
lambda = l_crit;
if ~isempty(DAE.V)
if CPF.show
fm_disp(['Continuation Power Flow completed in ', ...
num2str(toc),' s'],1);
end
maxgqerr = max(abs(DAE.gq));
maxgperr = max(abs(DAE.gp));
DAE.gq = zeros(Bus.n,1);
DAE.gp = zeros(Bus.n,1);
fm_call('pq');
glambda(Demand,lambda)
% load powers
Bus.Pl = DAE.gp;
Bus.Ql = DAE.gq;
% total bus injections
Vc = DAE.V.*exp(jay*DAE.a);
S = Vc.*conj(Line.Y*Vc);
% generated powers
Bus.Qg = imag(S) + Bus.Ql;
Bus.Pg = real(S) + Bus.Pl;
DAE.gq = maxgqerr*ones(Bus.n,1);
DAE.gp = maxgperr*ones(Bus.n,1);
if (Settings.showlf == 1 & CPF.show) | Fig.stat
SDbus = [Supply.bus;Demand.bus];
report = cell(2+length(dl_dp),1);
report{1,1} = ['Lambda_max = ', fvar(lambda_crit,9)];
report{2,1} = ' ';
for i = 1:length(dl_dp)
report{2+i,1} = ['d lambda / d P ', ...
Varname.bus{SDbus(i)},' = ', ...
fvar(dl_dp(i),9)];
end
fm_stat(report);
end
if CPF.show & Fig.plot
fm_plotfig
end
end
% Reset of SW, PV and PQ structures
PQ = restore(PQ);
PV = restore(PV);
SW = restore(SW);
if noDem, Demand = restore(Demand); end
if noSup, Supply = restore(Supply); end
if Syn.n
Syn.pm = Pm_syn;
end
if CPF.show & Fig.main
set(Fig.main,'Pointer','arrow');
Settings.xlabel = ['Loading Parameter ',char(92),'lambda (p.u.)'];
Settings.tf = 1.2*lambda_crit;
end
fm_status('cpf','close')
%warning on
Settings.pv2pq = PV2PQ;
CPF.lambda = lambda_crit;
CPF.kg = kg;
SNB.init = 0;
LIB.init = 0;
CPF.init = 1;
OPF.init = 0;
% --------------------------------
function s = signperm(P)
% --------------------------------
[i,j,p] = find(sparse(P));
idx = find(i ~= j);
q = P(idx,idx);
s = det(q);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -