基于半不变量法的随机潮流计算模型及方法
参考文献：《基于半不变量法的随机潮流计算模型及方法》
摘 要: 提出了一种线性化随机潮流的计算模型，在牛顿-拉夫逊法的基础上，利用半不变量法对随机变量进行卷积运算，并运用 Gram-Charlier 级数展开式计算随机变量的分布，得到节点电压概率密度函数． 通过节点电压概率分析为运行分析提供指导，以及对 IEEE-30节点系统进行计算，验证了所提模型和方法的有效性．
关键词: 电力系统; 随机潮流; 半不变量; 概率密度函数
1 引言
电力系统潮流计算是电力系统运行分析的基础,其方法包含确定性方法和不确定性潮流计算方法两大类.确定性潮流分析方法,即给定网络拓扑结构、元件的参数、节点负荷、发电机出力等参数，求解各节点电压及支路潮流的确定值。但由于实际运行中部分参数的不确定性，使得确定性潮流分析方法受到限制.为了准确模拟参数的多样性,也采用多场景计算方法,但这是对参数多种取值的枚举,计算量大,而且很难反映系统整体情况。因此,将参数的不确定性表示方法加入到潮流方程中,构造了不确定性潮流方程,形成了多种不确定性潮流分析方法,如随机潮流、盲数潮流、模糊潮流等。本文利用了随机潮流模型，改进了确定性潮流计算的不足，解决了潮流问题中的参数不确定性问题。
随机潮流计算运用概率统计方法处理系统运行中的随机变化因素,给出了系统运行电压、支路潮流的概率分布情况，更深刻地揭示了系统运行状况、系统中存在的问题、薄弱环节等信息，为规划与运行决策提供了更完整的信息 。因此,在电网规划中进行随机潮流计算显得尤为重要.此外在随机潮流计算中,本文利用了 Matlab 软件，大大缩短了程序编写和实现的过程,加快了程序的运算速度,提高了运算精度。
2 半不变量和 Cham-charlier 级数
半不变量是随机变量的一种数字特征，它可以由不高于相应阶次的随机变量的各阶矩求得。随机变量的各阶中心矩表示为:
可以得到各阶半不变量 Ki与中心矩 Mi的关系式为:
半不变量有以下重要性质:
由此可体现出半不变量的相加性，即独立随机变量之和的各阶半不变量等于各随机变量的各阶半不变量之和。
利用半不变量求和的运算代替卷积积分运算的一个关键问题是:已知矩量或半不变量,如何求出随机变量的概率密度函数.常用的方法是把概率密度函数展开成级数，用半不变量表示级数展开式的系数，本文采用 Gram-Charlier 级数展开。
3 计算步骤
由以上的理论分析可知，随机潮流计算的步骤可描述如图 1 所示。
4 算例
本文利用 Matlab 对 IEEE33 节点系统进行潮流计算分析.IEEE-33 节点系统的电网结构如图所示。
5 程序运行结果
1）节点6
2）节点14
3）节点22
4）节点26
6 程序
1）主函数
%% 基于半不变量法的随机潮流计算
clear
clc
close all
%% 半不变量法随机潮流---------------------------------------------CM-----------------------------------------------------
shuju=data_ieee30;
pv=find(shuju.bus(:,2)==2);
pq=find(shuju.bus(:,2)==1);
isb=find(shuju.bus(:,2)==3); %平衡节点
npv=length(pv);
npq=length(pq);
nb=size(shuju.bus,1); %节点数
mu_P_G=zeros(nb,1); %%节点发电机有功功率
mu_P_G(shuju.gen(:,1))=shuju.gen(:,2)/100;
mu_Q_G=zeros(nb,1); %%节点发电机无功功率
mu_Q_G(shuju.gen(:,1))=shuju.gen(:,3)/100;
mu_P_L=shuju.bus(:,3)/100; %%节点负荷有功功率
mu_Q_L=shuju.bus(:,4)/100; %%节点负荷无功功率
sigma_P_L=mu_P_L*0.3;
sigma_Q_L=mu_Q_L*0.3;
mu_P=mu_P_G-mu_P_L; %%节点注入有功
mu_Q=mu_Q_G-mu_Q_L; %%节点注入无功
%% -------计算节点注入功率的半不变量
jieshu=6;
gama_p=zeros(nb,jieshu); %注入有功的半不变量
gama_p(:,2)=sigma_P_L.^2;
gama_q=zeros(nb,jieshu); %注入无功的半不变量
gama_q(:,2)=sigma_Q_L.^2;
%-平衡节点及pv节点的Q置为0？
gama_p(isb,2)=0;
gama_q(isb,2)=0;
gama_q(pv,2)=0;
gama_pq=[gama_p;gama_q];
%% --------潮流计算形成2n阶雅克比矩阵
[basemva bus gen branch success et dSbus_dVm dSbus_dVa]=runpf(shuju);
H=real(dSbus_dVa);
M=imag(dSbus_dVa);
N=real(dSbus_dVm);
L=imag(dSbus_dVm);
Jacco2=[H N
M L]; %N=dP/dU L=dQ/dU
%-修正平衡节点和pv节点对应的系数
Jacco=Jacco2;
Jacco(isb,:)=0;
Jacco(:,isb)=0;
Jacco(isb,isb)=1;
Jacco(isb+nb,:)=0;
Jacco(:,isb+nb)=0;
Jacco(isb+nb,isb+nb)=1;
Jacco(pv+nb,:)=0;
Jacco(:,pv+nb)=0;
Jacco=Jacco+sparse(nb+pv, nb+pv, ones(npv,1), 2*nb, 2*nb);
%% ----------计算节点电压的半不变量
s1=Jacco\eye(size(Jacco));
gama_v=zeros(2*nb,jieshu); %电压相角、幅值
for i=1:jieshu
gama_v(:,i)=(s1.^i)*gama_pq(:,i);
end
gama_vm=gama_v(nb+pq,:); %pq节点电压幅值的半不变量
%-----------形成支路潮流的灵敏度矩阵
li=shuju.branch(:,1);
lj=shuju.branch(:,2);
b=size(shuju.branch,1); %线路数
G11=zeros(b,nb); %线路有功对角度的倒数
G12=zeros(b,nb); %线路有功对电压的导数
G21=zeros(b,nb);
G22=zeros(b,nb);
for i=1:b
G11(i,li(i))=-H(li(i),lj(i));
G11(i,lj(i))=H(li(i),lj(i));
G12(i,li(i))=0.01*2*branch(i,14)/bus(li(i),8)-N(li(i),lj(i)); %%化为标幺值
G12(i,lj(i))=N(li(i),lj(i));
G21(i,li(i))=-M(li(i),lj(i));
G21(i,lj(i))=M(li(i),lj(i));
G22(i,li(i))=0.01*2*branch(i,15)/bus(li(i),8)-H(li(i),lj(i)); %%化为标幺值
G22(i,lj(i))=H(li(i),lj(i));
end
G=[G11 G12;G21 G22];
t1=G*s1;
%------------计算支路潮流的半不变量
gama_xianlu=zeros(2*b,jieshu);
for i=1:jieshu
gama_xianlu(:,i)=(t1.^i)*gama_pq(:,i);
end
gama_xianlu_p=gama_xianlu(1:b,:);
gama_xianlu_q=gama_xianlu(b+1:end,:);
%% --------------用半不变量求Gram-charlie展开系数
g_vm=zeros(size(gama_vm));
g_xianlu_p=zeros(size(gama_xianlu_p));
g_xianlu_q=zeros(size(gama_xianlu_q));
for i=1:jieshu
g_vm(:,i)=gama_vm(:,i)./(gama_vm(:,2).^(i/2));
g_xianlu_p(:,i)=gama_xianlu_p(:,i)./(gama_xianlu_p(:,2).^(i/2));
g_xianlu_q(:,i)=gama_xianlu_q(:,i)./(gama_xianlu_q(:,2).^(i/2));
end
%---电压展开系数C
c_vm=zeros(size(g_vm));
c_vm(:,3)=-g_vm(:,3);
c_vm(:,4)=g_vm(:,4);
c_vm(:,5)=-g_vm(:,5);
c_vm(:,6)=g_vm(:,6)+10*g_vm(:,3).^2;
%---线路有功展开系数C
c_xianlu_p=zeros(size(g_xianlu_p));
c_xianlu_p(:,3)=-g_xianlu_p(:,3);
c_xianlu_p(:,4)=g_xianlu_p(:,4);
c_xianlu_p(:,5)=-g_xianlu_p(:,5);
c_xianlu_p(:,6)=g_xianlu_p(:,6)+10*g_xianlu_p(:,3).^2;
%---线路无功展开系数C
c_xianlu_q=zeros(size(g_xianlu_q));
c_xianlu_q(:,3)=-g_xianlu_q(:,3);
c_xianlu_q(:,4)=g_xianlu_q(:,4);
c_xianlu_q(:,5)=-g_xianlu_q(:,5);
c_xianlu_q(:,6)=g_xianlu_q(:,6)+10*g_xianlu_q(:,3).^2;
%% -pq节点电压幅值期望及方差、线路功率期望及方差
mu_vm=bus(pq,8);
sigma_vm=sqrt(gama_vm(:,2));
mu_xianlu_p=branch(:,14)/100;
sigma_xianlu_p=sqrt(gama_xianlu_p(:,2));
mu_xianlu_q=branch(:,15)/100;
sigma_xianlu_q=sqrt(gama_xianlu_q(:,2));
%pq节点电压概率密度
for i=1:npq
x=0.98:0.0001:1.05;
x1=(x-mu_vm(i))/sigma_vm(i);
pdf_vm(i,:)=c2pdf(x1,c_vm(i,:),sigma_vm(i));
end
figure
plot(x,pdf_vm(2,:),'r--')
%线路1有功
px=0.2:0.0001:0.9;
px1=(px-mu_xianlu_p(1))/sigma_xianlu_p(1); %线路1有功
pdf_xianlu1_p=c2pdf(px1,c_xianlu_p(1,:),sigma_xianlu_p(1));
figure
plot(px,pdf_xianlu1_p,'r--')
legend('半不变量法')
title('线路1有功概率密度曲线')
%------------绘制各节点的概率分布
for i=1:npq
figure
% ksdensity(v_mc(pq(i),:))
% hold on
plot(x,pdf_vm(i,:),'r--')
legend('半不变量法')
title(['节点',num2str(pq(i)),'概率密度曲线'])
end
% % %---------由展开系数得到各点对应的概率密度
function pdf_x=c2pdf(x,c,sigma)
pdf_x=normpdf(x).*(1+c(3)*(3*x-x.^3)/factorial(3)+c(4)*(x.^4-6*x.^2+3)/factorial(4)+...
c(5)*(-x.^5+10*x.^3-15*x)/factorial(5)+c(6)*(x.^6-15*x.^4+45*x.^2-15)/factorial(6))/sigma;
end
%----------由展开系数得到累计分布
function cdf_x=c2cdf(x,c)
cdf_x=normcdf(x)+normpdf(x).*(c(3)*(x.^2-1)/factorial(3)+c(4)*(-x.^3+3*x)/factorial(4)+...
c(5)*(x.^4-6*x.^2+3)/factorial(5)+c(6)*(-x.^5+10*x.^3-15*x)/factorial(6));
end
2）子函数
function Sd = makeSdzip(baseMVA, bus, mpopt)
%MAKESDZIP Builds vectors of nominal complex bus power demands for ZIP loads.
% SD = MAKESDZIP(BASEMVA, BUS, MPOPT) returns a struct with three fields,
% each an nb x 1 vectors. The fields 'z', 'i' and 'p' correspond to the
% nominal p.u. complex power (at 1 p.u. voltage magnitude) of the constant
% impedance, constant current, and constant power portions, respectively of
% the ZIP load model.
%
% Example:
% Sd = makeSdzip(baseMVA, bus, mpopt);
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
if nargin < 3
mpopt = [];
end
if ~isempty(mpopt) && ~isempty(mpopt.exp.sys_wide_zip_loads.pw)
if any(size(mpopt.exp.sys_wide_zip_loads.pw) ~= [1 3])
error('makeSdzip: ''exp.sys_wide_zip_loads.pw'' must be a 1 x 3 vector');
end
if abs(sum(mpopt.exp.sys_wide_zip_loads.pw) - 1) > eps
error('makeSdzip: elements of ''exp.sys_wide_zip_loads.pw'' must sum to 1');
end
pw = mpopt.exp.sys_wide_zip_loads.pw;
else
pw = [1 0 0];
end
if ~isempty(mpopt) && ~isempty(mpopt.exp.sys_wide_zip_loads.qw)
if any(size(mpopt.exp.sys_wide_zip_loads.qw) ~= [1 3])
error('makeSdzip: ''exp.sys_wide_zip_loads.qw'' must be a 1 x 3 vector');
end
if abs(sum(mpopt.exp.sys_wide_zip_loads.qw) - 1) > eps
error('makeSdzip: elements of ''exp.sys_wide_zip_loads.qw'' must sum to 1');
end
qw = mpopt.exp.sys_wide_zip_loads.qw;
else
qw = pw;
end
Sd.z = (bus(:, PD) * pw(3) + 1j * bus(:, QD) * qw(3)) / baseMVA;
Sd.i = (bus(:, PD) * pw(2) + 1j * bus(:, QD) * qw(2)) / baseMVA;
Sd.p = (bus(:, PD) * pw(1) + 1j * bus(:, QD) * qw(1)) / baseMVA;
function [Sbus, dSbus_dVm] = makeSbus(baseMVA, bus, gen, mpopt, Vm, Sg)
%MAKESBUS Builds the vector of complex bus power injections.
% SBUS = MAKESBUS(BASEMVA, BUS, GEN)
% SBUS = MAKESBUS(BASEMVA, BUS, GEN, MPOPT, VM)
% SBUS = MAKESBUS(BASEMVA, BUS, GEN, MPOPT, VM, SG)
% returns the vector of complex bus power injections, that is, generation
% minus load. Power is expressed in per unit. If the MPOPT and VM arguments
% are present it evaluates any ZIP loads based on the provided voltage
% magnitude vector. If VM is empty, it assumes nominal voltage. If SG is
% provided, it is a complex ng x 1 vector of generator power injections in
%% define named indices into bus, gen matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
%% default inputs
if nargin < 5
Vm = [];
if nargin < 4
mpopt = [];
end
end
nb = size(bus, 1);
%% get load parameters
Sd = makeSdzip(baseMVA, bus, mpopt);
if nargout == 2
Sbus = [];
if isempty(Vm)
dSbus_dVm = sparse(nb, nb);
else
dSbus_dVm = -(spdiags(Sd.i + 2 * Vm .* Sd.z, 0, nb, nb));
end
else
%% compute per-bus generation in p.u.
on = find(gen(:, GEN_STATUS) > 0); %% which generators are on?
gbus = gen(on, GEN_BUS); %% what buses are they at?
ngon = size(on, 1);
Cg = sparse(gbus, (1:ngon)', 1, nb, ngon); %% connection matrix
%% element i, j is 1 if
%% gen on(j) at bus i is ON
if nargin > 5 && ~isempty(Sg)
Sbusg = Cg * Sg(on);
else
Sbusg = Cg * (gen(on, PG) + 1j * gen(on, QG)) / baseMVA;
end
%% compute per-bus loads in p.u.
if isempty(Vm)
Vm = ones(nb, 1);
end
Sbusd = Sd.p + Sd.i .* Vm + Sd.z .* Vm.^2;
%% form net complex bus power injection vector
%% (power injected by generators + power injected by loads)
Sbus = Sbusg - Sbusd;
end
function [MVAbase, bus, gen, branch, success, et, dSbus_dVm, dSbus_dVa] = ...
runpf(casedata, mpopt, fname, solvedcase)
%RUNPF Runs a power flow.
% [RESULTS, SUCCESS] = RUNPF(CASEDATA, MPOPT, FNAME, SOLVEDCASE)
%
% Runs a power flow (full AC Newton's method by default), optionally
% returning a RESULTS struct and SUCCESS flag.
%
% Inputs (all are optional):
% CASEDATA : either a MATPOWER case struct or a string containing
% the name of the file with the case data (default is 'case9')
% (see also CASEFORMAT and LOADCASE)
% MPOPT : MATPOWER options struct to override default options
% can be used to specify the solution algorithm, output options
% termination tolerances, and more (see also MPOPTION).
% FNAME : name of a file to which the pretty-printed output will
% be appended
% SOLVEDCASE : name of file to which the solved case will be saved
% in MATPOWER case format (M-file will be assumed unless the
% specified name ends with '.mat')
%
% Outputs (all are optional):
% RESULTS : results struct, with the following fields:
% (all fields from the input MATPOWER case, i.e. bus, branch,
% gen, etc., but with solved voltages, power flows, etc.)
% order - info used in external <-> internal data conversion
% et - elapsed time in seconds
% success - success flag, 1 = succeeded, 0 = failed
% SUCCESS : the success flag can additionally be returned as
% a second output argument
%
% Calling syntax options:
% results = runpf;
% results = runpf(casedata);
% results = runpf(casedata, mpopt);
% results = runpf(casedata, mpopt, fname);
% results = runpf(casedata, mpopt, fname, solvedcase);
% [results, success] = runpf(...);
%
% Alternatively, for compatibility with previous versions of MATPOWER,
% some of the results can be returned as individual output arguments:
%
% [baseMVA, bus, gen, branch, success, et] = runpf(...);
%
% If the pf.enforce_q_lims option is set to true (default is false) then, if
% any generator reactive power limit is violated after running the AC power
% flow, the corresponding bus is converted to a PQ bus, with Qg at the
% limit, and the case is re-run. The voltage magnitude at the bus will
% deviate from the specified value in order to satisfy the reactive power
% limit. If the reference bus is converted to PQ, the first remaining PV
% bus will be used as the slack bus for the next iteration. This may
% result in the real power output at this generator being slightly off
% from the specified values.
%
% Examples:
% results = runpf('case30');
% results = runpf('case30', mpoption('pf.enforce_q_lims', 1));
%
% See also RUNDCPF.
% MATPOWER
% Copyright (c) 1996-2016 by Power System Engineering Research Center (PSERC)
% by Ray Zimmerman, PSERC Cornell
% Enforcing of generator Q limits inspired by contributions
% from Mu Lin, Lincoln University, New Zealand (1/14/05).
%
% This file is part of MATPOWER.
% Covered by the 3-clause BSD License (see LICENSE file for details).
% See http://www.pserc.cornell.edu/matpower/ for more info.
%%----- initialize -----
%% define named indices into bus, gen, branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
%% default arguments
if nargin < 4
solvedcase = ''; %% don't save solved case
if nargin < 3
fname = ''; %% don't print results to a file
if nargin < 2
mpopt = mpoption; %% use default options
if nargin < 1
casedata = 'case9'; %% default data file is 'case9.m'
end
end
end
end
%% options
qlim = mpopt.pf.enforce_q_lims; %% enforce Q limits on gens?
dc = strcmp(upper(mpopt.model), 'DC'); %% use DC formulation?
%% read data
mpc = loadcase(casedata);
%% add zero columns to branch for flows if needed
if size(mpc.branch,2) < QT
mpc.branch = [ mpc.branch zeros(size(mpc.branch, 1), QT-size(mpc.branch,2)) ];
end
%% convert to internal indexing
mpc = ext2int(mpc);
[baseMVA, bus, gen, branch] = deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch);
%% get bus index lists of each type of bus
[ref, pv, pq] = bustypes(bus, gen);
%% generator info
on = find(gen(:, GEN_STATUS) > 0); %% which generators are on?
gbus = gen(on, GEN_BUS); %% what buses are they at?
%%----- run the power flow -----
t0 = clock;
if mpopt.verbose > 0
v = mpver('all');
fprintf('\nMATPOWER Version %s, %s', v.Version, v.Date);
end
if dc %% DC formulation
if mpopt.verbose > 0
fprintf(' -- DC Power Flow\n');
end
%% initial state
Va0 = bus(:, VA) * (pi/180);
%% build B matrices and phase shift injections
[B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch);
%% compute complex bus power injections (generation - load)
%% adjusted for phase shifters and real shunts
Pbus = real(makeSbus(baseMVA, bus, gen)) - Pbusinj - bus(:, GS) / baseMVA;
%% "run" the power flow
[Va, success] = dcpf(B, Pbus, Va0, ref, pv, pq);
%% update data matrices with solution
branch(:, [QF, QT]) = zeros(size(branch, 1), 2);
branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;
branch(:, PT) = -branch(:, PF);
bus(:, VM) = ones(size(bus, 1), 1);
bus(:, VA) = Va * (180/pi);
%% update Pg for slack generator (1st gen at ref bus)
%% (note: other gens at ref bus are accounted for in Pbus)
%% Pg = Pinj + Pload + Gs
%% newPg = oldPg + newPinj - oldPinj
refgen = zeros(size(ref));
for k = 1:length(ref)
temp = find(gbus == ref(k));
refgen(k) = on(temp(1));
end
gen(refgen, PG) = gen(refgen, PG) + (B(ref, :) * Va - Pbus(ref)) * baseMVA;
else %% AC formulation
alg = upper(mpopt.pf.alg);
if mpopt.verbose > 0
switch alg
case 'NR'
solver = 'Newton';
case 'FDXB'
solver = 'fast-decoupled, XB';
case 'FDBX'
solver = 'fast-decoupled, BX';
case 'GS'
solver = 'Gauss-Seidel';
otherwise
solver = 'unknown';
end
fprintf(' -- AC Power Flow (%s)\n', solver);
end
%% initial state
% V0 = ones(size(bus, 1), 1); %% flat start
V0 = bus(:, VM) .* exp(sqrt(-1) * pi/180 * bus(:, VA));
vcb = ones(size(V0)); %% create mask of voltage-controlled buses
vcb(pq) = 0; %% exclude PQ buses
k = find(vcb(gbus)); %% in-service gens at v-c buses
V0(gbus(k)) = gen(on(k), VG) ./ abs(V0(gbus(k))).* V0(gbus(k));
if qlim
ref0 = ref; %% save index and angle of
Varef0 = bus(ref0, VA); %% original reference bus(es)
limited = []; %% list of indices of gens @ Q lims
fixedQg = zeros(size(gen, 1), 1); %% Qg of gens at Q limits
end
%% build admittance matrices
[Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
repeat = 1;
while (repeat)
%% function for computing V dependent complex bus power injections
%% (generation - load)
Sbus = @(Vm)makeSbus(baseMVA, bus, gen, mpopt, Vm);
%% run the power flow
switch alg
case 'NR'
[V, success, iterations, dSbus_dVm, dSbus_dVa] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt);
case {'FDXB', 'FDBX'}
[Bp, Bpp] = makeB(baseMVA, bus, branch, alg);
[V, success, iterations] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt);
case 'GS'
if (~isempty(mpopt.exp.sys_wide_zip_loads.pw) && ...
any(mpopt.exp.sys_wide_zip_loads.pw(2:3))) || ...
(~isempty(mpopt.exp.sys_wide_zip_loads.qw) && ...
any(mpopt.exp.sys_wide_zip_loads.qw(2:3)))
warning('runpf: Gauss-Seidel algorithm does not support ZIP load model. Converting to constant power loads.')
mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads', ...
struct('pw', [], 'qw', []));
end
[V, success, iterations] = gausspf(Ybus, Sbus([]), V0, ref, pv, pq, mpopt);
otherwise
error('runpf: Only Newton''s method, fast-decoupled, and Gauss-Seidel power flow algorithms currently implemented.');
end
%% update data matrices with solution
[bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt);
if qlim %% enforce generator Q limits
%% find gens with violated Q constraints
mx = find( gen(:, GEN_STATUS) > 0 ...
& gen(:, QG) > gen(:, QMAX) + mpopt.opf.violation );
mn = find( gen(:, GEN_STATUS) > 0 ...
& gen(:, QG) < gen(:, QMIN) - mpopt.opf.violation );
if ~isempty(mx) || ~isempty(mn) %% we have some Q limit violations
%% first check for INFEASIBILITY (all remaining gens violating)
infeas = union(mx', mn')'; %% transposes handle fact that
%% union of scalars is a row vector
remaining = find( gen(:, GEN_STATUS) > 0 & ...
( bus(gen(:, GEN_BUS), BUS_TYPE) == PV | ...
bus(gen(:, GEN_BUS), BUS_TYPE) == REF ));
if length(infeas) == length(remaining) && all(infeas == remaining)
if mpopt.verbose
fprintf('All %d remaining gens exceed their Q limits : INFEASIBLE PROBLEM\n', length(infeas));
end
success = 0;
break;
end
%% one at a time?
if qlim == 2 %% fix largest violation, ignore the rest
[junk, k] = max([gen(mx, QG) - gen(mx, QMAX);
gen(mn, QMIN) - gen(mn, QG)]);
if k > length(mx)
mn = mn(k-length(mx));
mx = [];
else
mx = mx(k);
mn = [];
end
end
if mpopt.verbose && ~isempty(mx)
fprintf('Gen %d at upper Q limit, converting to PQ bus\n', mx);
end
if mpopt.verbose && ~isempty(mn)
fprintf('Gen %d at lower Q limit, converting to PQ bus\n', mn);
end
%% save corresponding limit values
fixedQg(mx) = gen(mx, QMAX);
fixedQg(mn) = gen(mn, QMIN);
mx = [mx;mn];
%% convert to PQ bus
gen(mx, QG) = fixedQg(mx); %% set Qg to binding limit
gen(mx, GEN_STATUS) = 0; %% temporarily turn off gen,
for i = 1:length(mx) %% (one at a time, since
bi = gen(mx(i), GEN_BUS); %% they may be at same bus)
bus(bi, [PD,QD]) = ... %% adjust load accordingly,
bus(bi, [PD,QD]) - gen(mx(i), [PG,QG]);
end
if length(ref) > 1 && any(bus(gen(mx, GEN_BUS), BUS_TYPE) == REF)
error('runpf: Sorry, MATPOWER cannot enforce Q limits for slack buses in systems with multiple slacks.');
end
bus(gen(mx, GEN_BUS), BUS_TYPE) = PQ; %% & set bus type to PQ
%% update bus index lists of each type of bus
ref_temp = ref;
[ref, pv, pq] = bustypes(bus, gen);
%% previous line can modify lists to select new REF bus
%% if there was none, so we should update bus with these
%% just to keep them consistent
if ref ~= ref_temp
bus(ref, BUS_TYPE) = REF;
bus( pv, BUS_TYPE) = PV;
if mpopt.verbose
fprintf('Bus %d is new slack bus\n', ref);
end
end
limited = [limited; mx];
else
repeat = 0; %% no more generator Q limits violated
end
else
repeat = 0; %% don't enforce generator Q limits, once is enough
end
end
if qlim && ~isempty(limited)
%% restore injections from limited gens (those at Q limits)
gen(limited, QG) = fixedQg(limited); %% restore Qg value,
for i = 1:length(limited) %% (one at a time, since
bi = gen(limited(i), GEN_BUS); %% they may be at same bus)
bus(bi, [PD,QD]) = ... %% re-adjust load,
bus(bi, [PD,QD]) + gen(limited(i), [PG,QG]);
end
gen(limited, GEN_STATUS) = 1; %% and turn gen back on
if ref ~= ref0
%% adjust voltage angles to make original ref bus correct
bus(:, VA) = bus(:, VA) - bus(ref0, VA) + Varef0;
end
end
end
mpc.et = etime(clock, t0);
mpc.success = success;
%%----- output results -----
%% convert back to original bus numbering & print results
[mpc.bus, mpc.gen, mpc.branch] = deal(bus, gen, branch);
results = int2ext(mpc);
%% zero out result fields of out-of-service gens & branches
if ~isempty(results.order.gen.status.off)
results.gen(results.order.gen.status.off, [PG QG]) = 0;
end
if ~isempty(results.order.branch.status.off)
results.branch(results.order.branch.status.off, [PF QF PT QT]) = 0;
end
if fname
[fd, msg] = fopen(fname, 'at');
if fd == -1
error(msg);
else
if mpopt.out.all == 0
printpf(results, fd, mpoption(mpopt, 'out.all', -1));
else
printpf(results, fd, mpopt);
end
fclose(fd);
end
end
printpf(results, 1, mpopt);
%% save solved case
if solvedcase
savecase(solvedcase, results);
end
if nargout == 1 || nargout == 2
MVAbase = results;
bus = success;
elseif nargout > 2
[MVAbase, bus, gen, branch, et] = ...
deal(results.baseMVA, results.bus, results.gen, results.branch, results.et);
% else %% don't define MVAbase, so it doesn't print anything
end
function [V, converged, i, dSbus_dVm, dSbus_dVa] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt)
%NEWTONPF Solves the power flow using a full Newton's method.
% [V, CONVERGED, I] = NEWTONPF(YBUS, SBUS, V0, REF, PV, PQ, MPOPT)
% solves for bus voltages given the full system admittance matrix (for
% all buses), the complex bus power injection vector (for all buses),
% the initial vector of complex bus voltages, and column vectors with
% the lists of bus indices for the swing bus, PV buses, and PQ buses,
% respectively. The bus voltage vector contains the set point for
% generator (including ref bus) buses, and the reference angle of the
% swing bus, as well as an initial guess for remaining magnitudes and
% angles. MPOPT is a MATPOWER options struct which can be used to
% set the termination tolerance, maximum number of iterations, and
% output options (see MPOPTION for details). Uses default options if
% this parameter is not given. Returns the final complex voltages, a
% flag which indicates whether it converged or not, and the number of
% iterations performed.
%
% See also RUNPF.
% MATPOWER
% Copyright (c) 1996-2016 by Power System Engineering Research Center (PSERC)
% by Ray Zimmerman, PSERC Cornell
%
% This file is part of MATPOWER.
% Covered by the 3-clause BSD License (see LICENSE file for details).
% See http://www.pserc.cornell.edu/matpower/ for more info.
%% default arguments
if nargin < 7
mpopt = mpoption;
end
%% options
tol = mpopt.pf.tol;
max_it = mpopt.pf.nr.max_it;
%% initialize
converged = 0;
i = 0;
V = V0;
Va = angle(V);
Vm = abs(V);
%% set up indexing for updating V
npv = length(pv);
npq = length(pq);
j1 = 1; j2 = npv; %% j1:j2 - V angle of pv buses
j3 = j2 + 1; j4 = j2 + npq; %% j3:j4 - V angle of pq buses
j5 = j4 + 1; j6 = j4 + npq; %% j5:j6 - V mag of pq buses
%% evaluate F(x0)
mis = V .* conj(Ybus * V) - Sbus(Vm);
F = [ real(mis([pv; pq]));
imag(mis(pq)) ];
%% check tolerance
normF = norm(F, inf);
if mpopt.verbose > 1
fprintf('\n it max P & Q mismatch (p.u.)');
fprintf('\n---- ---------------------------');
fprintf('\n%3d %10.3e', i, normF);
end
if normF < tol
converged = 1;
if mpopt.verbose > 1
fprintf('\nConverged!\n');
end
end
%% do Newton iterations
while (~converged && i < max_it)
%% update iteration counter
i = i + 1;
%% evaluate Jacobian
[dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V);
[dummy, neg_dSd_dVm] = Sbus(Vm);
dSbus_dVm = dSbus_dVm - neg_dSd_dVm;
j11 = real(dSbus_dVa([pv; pq], [pv; pq]));
j12 = real(dSbus_dVm([pv; pq], pq));
j21 = imag(dSbus_dVa(pq, [pv; pq]));
j22 = imag(dSbus_dVm(pq, pq));
J = [ j11 j12;
j21 j22; ];
%% compute update step
dx = -(J \ F);
%% update voltage
if npv
Va(pv) = Va(pv) + dx(j1:j2);
end
if npq
Va(pq) = Va(pq) + dx(j3:j4);
Vm(pq) = Vm(pq) + dx(j5:j6);
end
V = Vm .* exp(1j * Va);
Vm = abs(V); %% update Vm and Va again in case
Va = angle(V); %% we wrapped around with a negative Vm
%% evalute F(x)
mis = V .* conj(Ybus * V) - Sbus(Vm);
F = [ real(mis(pv));
real(mis(pq));
imag(mis(pq)) ];
%% check for convergence
normF = norm(F, inf);
if mpopt.verbose > 1
fprintf('\n%3d %10.3e', i, normF);
end
if normF < tol
converged = 1;
if mpopt.verbose
fprintf('\nNewton''s method power flow converged in %d iterations.\n', i);
end
end
end
if mpopt.verbose
if ~converged
fprintf('\nNewton''s method power flow did not converge in %d iterations.\n', i);
end
end
function mpc = data_ieee30
%CASE_IEEE30 Power flow data for IEEE 30 bus test case.
% Please see CASEFORMAT for details on the case file format.
% This data was converted from IEEE Common Data Format
% (ieee30cdf.txt) on 15-Oct-2014 by cdf2matp, rev. 2393
% See end of file for warnings generated during conversion.
%
% Converted from IEEE CDF file from:
% http://www.ee.washington.edu/research/pstca/
%
% 08/20/93 UW ARCHIVE 100.0 1961 W IEEE 30 Bus Test Case
% MATPOWER
%% MATPOWER Case Format : Version 2
mpc.version = '2';
%%----- Power Flow Data -----%%
%% system MVA base
mpc.baseMVA = 100;
%% bus data
% bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin
mpc.bus = [
1 3 0 0 0 0 1 1.05 0 132 1 1.06 0.94;
2 2 21.7 12.7 0 0 1 1.045 -5.48 132 1 1.06 0.94;
3 1 2.4 1.2 0 0 1 1.021 -7.96 132 1 1.06 0.94;
4 1 7.6 1.6 0 0 1 1.012 -9.62 132 1 1.06 0.94;
5 2 94.2 19 0 0 1 1.01 -14.37 132 1 1.06 0.94;
6 1 0 0 0 0 1 1.01 -11.34 132 1 1.06 0.94;
7 1 22.8 10.9 0 0 1 1.002 -13.12 132 1 1.06 0.94;
8 2 30 30 0 0 1 1.01 -12.1 132 1 1.06 0.94;
9 1 0 0 0 0 1 1.051 -14.38 1 1 1.06 0.94;
10 1 5.8 2 0 19 1 1.045 -15.97 33 1 1.06 0.94;
11 2 0 0 0 0 1 1.05 -14.39 11 1 1.06 0.94;
12 1 11.2 7.5 0 0 1 1.057 -15.24 33 1 1.06 0.94;
13 2 0 0 0 0 1 1.05 -15.24 11 1 1.06 0.94;
14 1 6.2 1.6 0 0 1 1.042 -16.13 33 1 1.06 0.94;
15 1 8.2 2.5 0 0 1 1.038 -16.22 33 1 1.06 0.94;
16 1 3.5 1.8 0 0 1 1.045 -15.83 33 1 1.06 0.94;
17 1 9 5.8 0 0 1 1.04 -16.14 33 1 1.06 0.94;
18 1 3.2 0.9 0 0 1 1.028 -16.82 33 1 1.06 0.94;
19 1 9.5 3.4 0 0 1 1.026 -17 33 1 1.06 0.94;
20 1 2.2 0.7 0 0 1 1.03 -16.8 33 1 1.06 0.94;
21 1 17.5 11.2 0 0 1 1.033 -16.42 33 1 1.06 0.94;
22 1 0 0 0 0 1 1.033 -16.41 33 1 1.06 0.94;
23 1 3.2 1.6 0 0 1 1.027 -16.61 33 1 1.06 0.94;
24 1 8.7 6.7 0 4.3 1 1.021 -16.78 33 1 1.06 0.94;
25 1 0 0 0 0 1 1.017 -16.35 33 1 1.06 0.94;
26 1 3.5 2.3 0 0 1 1 -16.77 33 1 1.06 0.94;
27 1 0 0 0 0 1 1.023 -15.82 33 1 1.06 0.94;
28 1 0 0 0 0 1 1.007 -11.97 132 1 1.06 0.94;
29 1 2.4 0.9 0 0 1 1.003 -17.06 33 1 1.06 0.94;
30 1 10.6 1.9 0 0 1 0.992 -17.94 33 1 1.06 0.94;
];
%% generator data
% bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10 ramp_30 ramp_q apf
mpc.gen = [
1 260.2 -16.1 150 -20 1.05 100 1 200 50 0 0 0 0 0 0 0 0 0 0 0;
2 80 50 60 -20 1.045 100 1 80 20 0 0 0 0 0 0 0 0 0 0 0;
5 50 37 62.45 -15 1.01 100 1 50 15 0 0 0 0 0 0 0 0 0 0 0;
8 20 37.3 48.73 -15 1.01 100 1 35 10 0 0 0 0 0 0 0 0 0 0 0;
11 20 16.2 40 -10 1.05 100 1 30 10 0 0 0 0 0 0 0 0 0 0 0;
13 20 10.6 44.72 -15 1.05 100 1 40 12 0 0 0 0 0 0 0 0 0 0 0;
];
%% branch data
% fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax
mpc.branch = [
1 2 0.0192 0.0575 0.0528 0 0 0 0 0 1 -360 360;
1 3 0.0452 0.1652 0.0408 0 0 0 0 0 1 -360 360;
2 4 0.057 0.1737 0.0368 0 0 0 0 0 1 -360 360;
3 4 0.0132 0.0379 0.0084 0 0 0 0 0 1 -360 360;
2 5 0.0472 0.1983 0.0418 0 0 0 0 0 1 -360 360;
2 6 0.0581 0.1763 0.0374 0 0 0 0 0 1 -360 360;
4 6 0.0119 0.0414 0.009 0 0 0 0 0 1 -360 360;
5 7 0.046 0.116 0.0204 0 0 0 0 0 1 -360 360;
6 7 0.0267 0.082 0.017 0 0 0 0 0 1 -360 360;
6 8 0.012 0.042 0.009 0 0 0 0 0 1 -360 360;
6 9 0 0.208 0 0 0 0 0.978 0 1 -360 360;
6 10 0 0.556 0 0 0 0 0.969 0 1 -360 360;
9 11 0 0.208 0 0 0 0 0 0 1 -360 360;
9 10 0 0.11 0 0 0 0 0 0 1 -360 360;
4 12 0 0.256 0 0 0 0 0.932 0 1 -360 360;
12 13 0 0.14 0 0 0 0 0 0 1 -360 360;
12 14 0.1231 0.2559 0 0 0 0 0 0 1 -360 360;
12 15 0.0662 0.1304 0 0 0 0 0 0 1 -360 360;
12 16 0.0945 0.1987 0 0 0 0 0 0 1 -360 360;
14 15 0.221 0.1997 0 0 0 0 0 0 1 -360 360;
16 17 0.0524 0.1923 0 0 0 0 0 0 1 -360 360;
15 18 0.1073 0.2185 0 0 0 0 0 0 1 -360 360;
18 19 0.0639 0.1292 0 0 0 0 0 0 1 -360 360;
19 20 0.034 0.068 0 0 0 0 0 0 1 -360 360;
10 20 0.0936 0.209 0 0 0 0 0 0 1 -360 360;
10 17 0.0324 0.0845 0 0 0 0 0 0 1 -360 360;
10 21 0.0348 0.0749 0 0 0 0 0 0 1 -360 360;
10 22 0.0727 0.1499 0 0 0 0 0 0 1 -360 360;
21 22 0.0116 0.0236 0 0 0 0 0 0 1 -360 360;
15 23 0.1 0.202 0 0 0 0 0 0 1 -360 360;
22 24 0.115 0.179 0 0 0 0 0 0 1 -360 360;
23 24 0.132 0.27 0 0 0 0 0 0 1 -360 360;
24 25 0.1885 0.3292 0 0 0 0 0 0 1 -360 360;
25 26 0.2544 0.38 0 0 0 0 0 0 1 -360 360;
25 27 0.1093 0.2087 0 0 0 0 0 0 1 -360 360;
28 27 0 0.396 0 0 0 0 0.968 0 1 -360 360;
27 29 0.2198 0.4153 0 0 0 0 0 0 1 -360 360;
27 30 0.3202 0.6027 0 0 0 0 0 0 1 -360 360;
29 30 0.2399 0.4533 0 0 0 0 0 0 1 -360 360;
8 28 0.0636 0.2 0.0428 0 0 0 0 0 1 -360 360;
6 28 0.0169 0.0599 0.013 0 0 0 0 0 1 -360 360;
];
%%----- OPF Data -----%%
%% generator cost data
% 1 startup shutdown n x1 y1 ... xn yn
% 2 startup shutdown n c(n-1) ... c0
mpc.gencost = [
2 0 0 3 0.0384319754 20 0;
2 0 0 3 0.25 20 0;
2 0 0 3 0.01 40 0;
2 0 0 3 0.01 40 0;
2 0 0 3 0.01 40 0;
2 0 0 3 0.01 40 0;
];
%% bus names
mpc.bus_name = {
'Glen Lyn 132';
'Claytor 132';
'Kumis 132';
'Hancock 132';
'Fieldale 132';
'Roanoke 132';
'Blaine 132';
'Reusens 132';
'Roanoke 1.0';
'Roanoke 33';
'Roanoke 11';
'Hancock 33';
'Hancock 11';
'Bus 14 33';
'Bus 15 33';
'Bus 16 33';
'Bus 17 33';
'Bus 18 33';
'Bus 19 33';
'Bus 20 33';
'Bus 21 33';
'Bus 22 33';
'Bus 23 33';
'Bus 24 33';
'Bus 25 33';
'Bus 26 33';
'Cloverdle 33';
'Cloverdle132';
'Bus 29 33';
'Bus 30 33';
};
% Warnings from cdf2matp conversion:
%
% ***** check the title format in the first line of the cdf file.
% ***** Qmax = Qmin at generator at bus 1 (Qmax set to Qmin + 10)
% ***** MVA limit of branch 1 - 2 not given, set to 0
% ***** MVA limit of branch 1 - 3 not given, set to 0
% ***** MVA limit of branch 2 - 4 not given, set to 0
% ***** MVA limit of branch 3 - 4 not given, set to 0
% ***** MVA limit of branch 2 - 5 not given, set to 0
% ***** MVA limit of branch 2 - 6 not given, set to 0
% ***** MVA limit of branch 4 - 6 not given, set to 0
% ***** MVA limit of branch 5 - 7 not given, set to 0
% ***** MVA limit of branch 6 - 7 not given, set to 0
% ***** MVA limit of branch 6 - 8 not given, set to 0
% ***** MVA limit of branch 6 - 9 not given, set to 0
% ***** MVA limit of branch 6 - 10 not given, set to 0
% ***** MVA limit of branch 9 - 11 not given, set to 0
% ***** MVA limit of branch 9 - 10 not given, set to 0
% ***** MVA limit of branch 4 - 12 not given, set to 0
% ***** MVA limit of branch 12 - 13 not given, set to 0
% ***** MVA limit of branch 12 - 14 not given, set to 0
% ***** MVA limit of branch 12 - 15 not given, set to 0
% ***** MVA limit of branch 12 - 16 not given, set to 0
% ***** MVA limit of branch 14 - 15 not given, set to 0
% ***** MVA limit of branch 16 - 17 not given, set to 0
% ***** MVA limit of branch 15 - 18 not given, set to 0
% ***** MVA limit of branch 18 - 19 not given, set to 0
% ***** MVA limit of branch 19 - 20 not given, set to 0
% ***** MVA limit of branch 10 - 20 not given, set to 0
% ***** MVA limit of branch 10 - 17 not given, set to 0
% ***** MVA limit of branch 10 - 21 not given, set to 0
% ***** MVA limit of branch 10 - 22 not given, set to 0
% ***** MVA limit of branch 21 - 22 not given, set to 0
% ***** MVA limit of branch 15 - 23 not given, set to 0
% ***** MVA limit of branch 22 - 24 not given, set to 0
% ***** MVA limit of branch 23 - 24 not given, set to 0
% ***** MVA limit of branch 24 - 25 not given, set to 0
% ***** MVA limit of branch 25 - 26 not given, set to 0
% ***** MVA limit of branch 25 - 27 not given, set to 0
% ***** MVA limit of branch 28 - 27 not given, set to 0
% ***** MVA limit of branch 27 - 29 not given, set to 0
% ***** MVA limit of branch 27 - 30 not given, set to 0
% ***** MVA limit of branch 29 - 30 not given, set to 0
% ***** MVA limit of branch 8 - 28 not given, set to 0
% ***** MVA limit of branch 6 - 28 not given, set to 0
function opt = mpoption(varargin)
%MPOPTION Used to set and retrieve a MATPOWER options struct.
%
% OPT = MPOPTION
% Returns the default options struct.
%
% OPT = MPOPTION(OVERRIDES)
% Returns the default options struct, with some fields overridden
% by values from OVERRIDES, which can be a struct or the name of
% a function that returns a struct.
%T0 is an options struct it does nothing.
%
% OPT = MPOPTION(OPT0, OVERRIDES)
% Applies overrides to an existing set of options, OPT0, which
% can be an old-style options vector or an options struct.
%
% OPT = MPOPTION(OPT0, NAME1, VALUE1, NAME2, VALUE2, ...)
% Same as above except it uses the old-style options vector OPT0
% as a base instead of the old default options vector.
%
% OPT_VECTOR = MPOPTION(OPT, [])
% Creates and returns an old-style options vector from an
% options struct OPT.
%
% Note: The use of old-style MATPOWER options vectors and their
% names and values has been deprecated and will be removed
% in a future version of MATPOWER. Until then, all uppercase
% option names are not permitted for new top-level options.
%
% Examples:
% mpopt = mpoption('pf.alg', 'FDXB', 'pf.tol', 1e-4);
% mpopt = mpoption(mpopt, 'opf.dc.solver', 'CPLEX', 'verbose', 2);
%
%The currently defined options are as follows:
%
% name default description [options]
%---------------------- --------- ----------------------------------
%Model options:
% model 'AC' AC vs. DC power flow model
% [ 'AC' - use nonlinear AC model & corresponding algorithms/options ]
% [ 'DC' - use linear DC model & corresponding algorithms/options ]
%
%Power Flow options:
% pf.alg 'NR' AC power flow algorithm
% [ 'NR' - Newton's method ]
% [ 'FDXB' - Fast-Decoupled (XB version) ]
% [ 'FDBX' - Fast-Decoupled (BX version) ]
% [ 'GS' - Gauss-Seidel ]
% pf.tol 1e-8 termination tolerance on per unit
% P & Q mismatch
% pf.nr.max_it 10 maximum number of iterations for
% Newton's method
% pf.fd.max_it 30 maximum number of iterations for
% fast decoupled method
% pf.gs.max_it 1000 maximum number of iterations for
% Gauss-Seidel method
% pf.enforce_q_lims 0 enforce gen reactive power limits at
% expense of |V|
% [ 0 - do NOT enforce limits ]
% [ 1 - enforce limits, simultaneous bus type conversion ]
% [ 2 - enforce limits, one-at-a-time bus type conversion ]
%
%Continuation Power Flow options:
% cpf.parameterization 3 choice of parameterization
% [ 1 - natural ]
% [ 2 - arc length ]
% [ 3 - pseudo arc length ]
% cpf.stop_at 'NOSE' determins stopping criterion
% [ 'NOSE' - stop when nose point is reached ]
% [ 'FULL' - trace full nose curve ]
% [ <lam_stop> - stop upon reaching specified target lambda value ]
% cpf.step 0.05 continuation power flow step size
% cpf.adapt_step 0 toggle adaptive step size feature
% [ 0 - adaptive step size disabled ]
% [ 1 - adaptive step size enabled ]
% cpf.error_tol 1e-3 tolerance for adaptive step control
% cpf.step_min 1e-4 minimum allowed step size
% cpf.step_max 0.2 maximum allowed step size
% cpf.plot.level 0 control plotting of noze curve
% [ 0 - do not plot nose curve ]
% [ 1 - plot when completed ]
% [ 2 - plot incrementally at each iteration ]
% [ 3 - same as 2, with 'pause' at each iteration ]
% cpf.plot.bus <empty> index of bus whose voltage is to be
% plotted
% cpf.user_callback <empty> string or cell array of strings
% with names of user callback functions
% see 'help cpf_default_callback'
% cpf.user_callback_args <empty> struct passed to user-defined
% callback functions
%
%Optimal Power Flow options:
% name default description [options]
%---------------------- --------- ----------------------------------
% opf.ac.solver 'DEFAULT' AC optimal power flow solver
% [ 'DEFAULT' - choose solver based on availability in the following ]
% [ order: 'PDIPM', 'MIPS' ]
% [ 'MIPS' - MIPS, Matlab Interior Point Solver, primal/dual ]
% [ interior point method (pure Matlab) ]
% [ 'FMINCON' - MATLAB Optimization Toolbox, FMINCON ]
% [ 'IPOPT' - IPOPT, requires MEX interface to IPOPT solver ]
% [ available from: ]
% [ http://www.coin-or.org/projects/Ipopt.xml ]
% [ 'KNITRO' - KNITRO, requires MATLAB Optimization Toolbox and ]
% [ KNITRO libraries available from: http://www.ziena.com/]
% [ 'MINOPF' - MINOPF, MINOS-based solver, requires optional ]
% [ MEX-based MINOPF package, available from: ]
% [ http://www.pserc.cornell.edu/minopf/ ]
% [ 'PDIPM' - PDIPM, primal/dual interior point method, requires ]
% [ optional MEX-based TSPOPF package, available from: ]
% [ http://www.pserc.cornell.edu/tspopf/ ]
% [ 'SDPOPF' - SDPOPF, solver based on semidefinite relaxation of ]
% [ OPF problem, requires optional packages: ]
% [ SDP_PF, available in extras/sdp_pf ]
% [ YALMIP, available from: ]
% [ http://users.isy.liu.se/johanl/yalmip/ ]
% [ SDP solver such as SeDuMi, available from: ]
% [ http://sedumi.ie.lehigh.edu/ ]
% [ 'TRALM' - TRALM, trust region based augmented Langrangian ]
% [ method, requires TSPOPF (see 'PDIPM') ]
% opf.dc.solver 'DEFAULT' DC optimal power flow solver
% [ 'DEFAULT' - choose solver based on availability in the following ]
% [ order: 'GUROBI', 'CPLEX', 'MOSEK', 'OT', ]
% [ 'GLPK' (linear costs only), 'BPMPD', 'MIPS' ]
% [ 'MIPS' - MIPS, Matlab Interior Point Solver, primal/dual ]
% [ interior point method (pure Matlab) ]
% [ 'BPMPD' - BPMPD, requires optional MEX-based BPMPD_MEX package ]
% [ available from: http://www.pserc.cornell.edu/bpmpd/ ]
% [ 'CLP' - CLP, requires interface to COIN-OP LP solver ]
% [ available from:http://www.coin-or.org/projects/Clp.xml]
% [ 'CPLEX' - CPLEX, requires CPLEX solver available from: ]
% [ http://www.ibm.com/software/integration/ ... ]
% [ ... optimization/cplex-optimizer/ ]
% [ 'GLPK' - GLPK, requires interface to GLPK solver ]
% [ available from: http://www.gnu.org/software/glpk/ ]
% [ (GLPK does not work with quadratic cost functions) ]
% [ 'GUROBI' - GUROBI, requires Gurobi optimizer (v. 5+) ]
% [ available from: http://www.gurobi.com/ ]
% [ 'IPOPT' - IPOPT, requires MEX interface to IPOPT solver ]
% [ available from: ]
% [ http://www.coin-or.org/projects/Ipopt.xml ]
% [ 'MOSEK' - MOSEK, requires Matlab interface to MOSEK solver ]
% [ available from: http://www.mosek.com/ ]
% [ 'OT' - MATLAB Optimization Toolbox, QUADPROG, LINPROG ]
% opf.violation 5e-6 constraint violation tolerance
% opf.flow_lim 'S' quantity limited by branch flow
% constraints
% [ 'S' - apparent power flow (limit in MVA) ]
% [ 'P' - active power flow (limit in MW) ]
% [ 'I' - current magnitude (limit in MVA at 1 p.u. voltage) ]
% opf.ignore_angle_lim 0 angle diff limits for branches
% [ 0 - include angle difference limits, if specified ]
% [ 1 - ignore angle difference limits even if specified ]
% opf.init_from_mpc -1 specify whether to use current state
% in MATPOWER case to initialize OPF
% (currently supported only for Ipopt,
% Knitro and MIPS solvers)
% [ -1 - MATPOWER decides, based on solver/algorithm ]
% [ 0 - ignore current state when initializing OPF ]
% [ 1 - use current state to initialize OPF ]
% opf.return_raw_der 0 for AC OPF, return constraint and
% derivative info in results.raw
% (in fields g, dg, df, d2f) [ 0 or 1 ]
%
%Output options:
% name default description [options]
%---------------------- --------- ----------------------------------
% verbose 1 amount of progress info printed
% [ 0 - print no progress info ]
% [ 1 - print a little progress info ]
% [ 2 - print a lot of progress info ]
% [ 3 - print all progress info ]
% out.all -1 controls pretty-printing of results
% [ -1 - individual flags control what prints ]
% [ 0 - do not print anything (overrides individual flags, ignored ]
% [ for files specified as FNAME arg to runpf(), runopf(), etc.)]
% [ 1 - print everything (overrides individual flags) ]
% out.sys_sum 1 print system summary [ 0 or 1 ]
% out.area_sum 0 print area summaries [ 0 or 1 ]
% out.bus 1 print bus detail [ 0 or 1 ]
% out.branch 1 print branch detail [ 0 or 1 ]
% out.gen 0 print generator detail [ 0 or 1 ]
% out.lim.all -1 controls constraint info output
% [ -1 - individual flags control what constraint info prints ]
% [ 0 - no constraint info (overrides individual flags) ]
% [ 1 - binding constraint info (overrides individual flags) ]
% [ 2 - all constraint info (overrides individual flags) ]
% out.lim.v 1 control voltage limit info
% [ 0 - do not print ]
% [ 1 - print binding constraints only ]
% [ 2 - print all constraints ]
% [ (same options for OUT_LINE_LIM, OUT_PG_LIM, OUT_QG_LIM) ]
% out.lim.line 1 control line flow limit info
% out.lim.pg 1 control gen active power limit info
% out.lim.qg 1 control gen reactive pwr limit info
% out.force 0 print results even if success
% flag = 0 [ 0 or 1 ]
% out.suppress_detail -1 suppress all output but system summary
% [ -1 - suppress details for large systems (> 500 buses) ]
% [ 0 - do not suppress any output specified by other flags ]
% [ 1 - suppress all output except system summary section ]
% [ (overrides individual flags, but not out.all = 1) ]
%
%Solver specific options:
% name default description [options]
% ----------------------- --------- ----------------------------------
% MIPS:
% mips.linsolver '' linear system solver
% [ '' or '\' build-in backslash \ operator (e.g. x = A \ b) ]
% [ 'PARDISO' PARDISO solver (if available) ]
% mips.feastol 0 feasibility (equality) tolerance
% (set to opf.violation by default)
% mips.gradtol 1e-6 gradient tolerance
% mips.comptol 1e-6 complementary condition
% (inequality) tolerance
% mips.costtol 1e-6 optimality tolerance
% mips.max_it 150 maximum number of iterations
% mips.step_control 0 enable step-size cntrl [ 0 or 1 ]
% mips.sc.red_it 20 maximum number of reductions per
% iteration with step control
% mips.xi 0.99995 constant used in alpha updates*
% mips.sigma 0.1 centering parameter*
% mips.z0 1 used to initialize slack variables*
% mips.alpha_min 1e-8 returns "Numerically Failed" if
% either alpha parameter becomes
% smaller than this value*
% mips.rho_min 0.95 lower bound on rho_t*
% mips.rho_max 1.05 upper bound on rho_t*
% mips.mu_threshold 1e-5 KT multipliers smaller than this
% value for non-binding constraints
% are forced to zero
% mips.max_stepsize 1e10 returns "Numerically Failed" if the
% 2-norm of the reduced Newton step
% exceeds this value*
% * See the corresponding Appendix in the manual for details.
%
% CPLEX:
% cplex.lpmethod 0 solution algorithm for LP problems
% [ 0 - automatic: let CPLEX choose ]
% [ 1 - primal simplex ]
% [ 2 - dual simplex ]
% [ 3 - network simplex ]
% [ 4 - barrier ]
% [ 5 - sifting ]
% [ 6 - concurrent (dual, barrier, and primal) ]
% cplex.qpmethod 0 solution algorithm for QP problems
% [ 0 - automatic: let CPLEX choose ]
% [ 1 - primal simplex optimizer ]
% [ 2 - dual simplex optimizer ]
% [ 3 - network optimizer ]
% [ 4 - barrier optimizer ]
% cplex.opts <empty> see CPLEX_OPTIONS for details
% cplex.opt_fname <empty> see CPLEX_OPTIONS for details
% cplex.opt 0 see CPLEX_OPTIONS for details
%
% FMINCON:
% fmincon.alg 4 algorithm used by fmincon() for OPF
% for Opt Toolbox 4 and later
% [ 1 - active-set (not suitable for large problems) ]
% [ 2 - interior-point, w/default 'bfgs' Hessian approx ]
% [ 3 - interior-point, w/ 'lbfgs' Hessian approx ]
% [ 4 - interior-point, w/exact user-supplied Hessian ]
% [ 5 - interior-point, w/Hessian via finite differences ]
% [ 6 - sqp (not suitable for large problems) ]
% fmincon.tol_x 1e-4 termination tol on x
% fmincon.tol_f 1e-4 termination tol on f
% fmincon.max_it 0 maximum number of iterations
% [ 0 => default ]
%
% GUROBI:
% gurobi.method 0 solution algorithm (Method)
% [ -1 - automatic, let Gurobi decide ]
% [ 0 - primal simplex ]
% [ 1 - dual simplex ]
% [ 2 - barrier ]
% [ 3 - concurrent (LP only) ]
% [ 4 - deterministic concurrent (LP only) ]
% gurobi.timelimit Inf maximum time allowed (TimeLimit)
% gurobi.threads 0 max number of threads (Threads)
% gurobi.opts <empty> see GUROBI_OPTIONS for details
% gurobi.opt_fname <empty> see GUROBI_OPTIONS for details
% gurobi.opt 0 see GUROBI_OPTIONS for details
%
% IPOPT:
% ipopt.opts <empty> see IPOPT_OPTIONS for details
% ipopt.opt_fname <empty> see IPOPT_OPTIONS for details
% ipopt.opt 0 see IPOPT_OPTIONS for details
%
% KNITRO:
% knitro.tol_x 1e-4 termination tol on x
% knitro.tol_f 1e-4 termination tol on f
% knitro.opt_fname <empty> name of user-supplied native
% KNITRO options file that overrides
% all other options
% knitro.opt 0 if knitro.opt_fname is empty and
% knitro.opt is a non-zero integer N
% then knitro.opt_fname is auto-
% generated as:
% 'knitro_user_options_N.txt'
%
% LINPROG:
% linprog <empty> LINPROG options passed to
% OPTIMOPTIONS or OPTIMSET.
% see LINPROG in the Optimization
% Toolbox for details
%
% MINOPF:
% minopf.feastol 0 (1e-3) primal feasibility tolerance
% (set to opf.violation by default)
% minopf.rowtol 0 (1e-3) row tolerance
% minopf.xtol 0 (1e-4) x tolerance
% minopf.majdamp 0 (0.5) major damping parameter
% minopf.mindamp 0 (2.0) minor damping parameter
% minopf.penalty 0 (1.0) penalty parameter
% minopf.major_it 0 (200) major iterations
% minopf.minor_it 0 (2500) minor iterations
% minopf.max_it 0 (2500) iterations limit
% minopf.verbosity -1 amount of progress info printed
% [ -1 - controlled by 'verbose' option ]
% [ 0 - print nothing ]
% [ 1 - print only termination status message ]
% [ 2 - print termination status and screen progress ]
% [ 3 - print screen progress, report file (usually fort.9) ]
% minopf.core 0 (1200*nb + 2*(nb+ng)^2) memory allocation
% minopf.supbasic_lim 0 (2*nb + 2*ng) superbasics limit
% minopf.mult_price 0 (30) multiple price
%
% MOSEK:
% mosek.lp_alg 0 solution algorithm for LP problems
% (MSK_IPAR_OPTIMIZER)
% for MOSEK 7.x ... (see MOSEK_SYMBCON for a "better way")
% [ 0 - automatic: let MOSEK choose ]
% [ 1 - interior point ]
% [ 3 - primal simplex ]
% [ 4 - dual simplex ]
% [ 5 - primal dual simplex ]
% [ 6 - automatic simplex (MOSEK chooses which simplex method) ]
% [ 7 - network primal simplex ]
% [ 10 - concurrent ]
% mosek.max_it 0 (400) interior point max iterations
% (MSK_IPAR_INTPNT_MAX_ITERATIONS)
% mosek.gap_tol 0 (1e-8) interior point relative gap tol
% (MSK_DPAR_INTPNT_TOL_REL_GAP)
% mosek.max_time 0 (-1) maximum time allowed
% (MSK_DPAR_OPTIMIZER_MAX_TIME)
% mosek.num_threads 0 (1) max number of threads
% (MSK_IPAR_INTPNT_NUM_THREADS)
% mosek.opts <empty> see MOSEK_OPTIONS for details
% mosek.opt_fname <empty> see MOSEK_OPTIONS for details
% mosek.opt 0 see MOSEK_OPTIONS for details
%
% QUADPROG:
% quadprog <empty> QUADPROG options passed to
% OPTIMOPTIONS or OPTIMSET.
% see QUADPROG in the Optimization
% Toolbox for details
%
% TSPOPF:
% pdipm.feastol 0 feasibility (equality) tolerance
% (set to opf.violation by default)
% pdipm.gradtol 1e-6 gradient tolerance
% pdipm.comptol 1e-6 complementary condition
% (inequality) tolerance
% pdipm.costtol 1e-6 optimality tolerance
% pdipm.max_it 150 maximum number of iterations
% pdipm.step_control 0 enable step-size cntrl [ 0 or 1 ]
% pdipm.sc.red_it 20 maximum number of reductions per
% iteration with step control
% pdipm.sc.smooth_ratio 0.04 piecewise linear curve smoothing
% ratio
%
% tralm.feastol 0 feasibility tolerance
% (set to opf.violation by default)
% tralm.primaltol 5e-4 primal variable tolerance
% tralm.dualtol 5e-4 dual variable tolerance
% tralm.costtol 1e-5 optimality tolerance
% tralm.major_it 40 maximum number of major iterations
% tralm.minor_it 40 maximum number of minor iterations
% tralm.smooth_ratio 0.04 piecewise linear curve smoothing
% ratio
%
%Experimental Options:
% exp.sys_wide_zip_loads.pw <empty> 1 x 3 vector of active load fraction
% to be modeled as constant power,
% constant current and constant
% impedance, respectively, where
% <empty> means use [1 0 0]
% exp.sys_wide_zip_loads.qw <empty> same for reactive power, where
% <empty> means use same value as
% for 'pw'
% MATPOWER
% Copyright (c) 2013-2016 by Power System Engineering Research Center (PSERC)
% by Ray Zimmerman, PSERC Cornell
%
% This file is part of MATPOWER.
% Covered by the 3-clause BSD License (see LICENSE file for details).
% See http://www.pserc.cornell.edu/matpower/ for more info.
%% some constants
N = 124; %% number of options in old-style vector
v = mpoption_version; %% version number of MATPOWER options struct
%% initialize flags and arg counter
have_opt0 = 0; %% existing options struct or vector provided?
have_old_style_ov = 0; %% override options using old-style names?
return_old_style = 0; %% return value as old-style vector?
k = 1;
if nargin > 0
opt0 = varargin{k};
if (isstruct(opt0) && isfield(opt0, 'v')) || ...
(isnumeric(opt0) && size(opt0, 1) == N && size(opt0, 2) == 1)
have_opt0 = 1;
k = k + 1;
end
end
%% create base options vector to which overrides are made
if have_opt0
if isstruct(opt0) %% it's already a valid options struct
if DEBUG, fprintf('OPT0 is a valid options struct\n'); end
if opt0.v < v
%% convert older version to current version
opt_d = mpoption_default();
if opt0.v == 1 %% convert version 1 to 2
if isfield(opt_d, 'linprog')
opt0.lingprog = opt_d.linprog;
end
if isfield(opt_d, 'quadprog')
opt0.quadprog = opt_d.quadprog;
end
end
if opt0.v <= 2 %% convert version 2 to 3
opt0.out.suppress_detail = opt_d.out.suppress_detail;
end
%if opt0.v <= 3 %% convert version 3 to 4
%% new mips options were all optional, no conversion needed
%end
if opt0.v <= 4 %% convert version 4 to 5
opt0.opf.init_from_mpc = opt_d.opf.init_from_mpc;
end
if opt0.v <= 5 %% convert version 5 to 6
if isfield(opt_d, 'clp')
opt0.clp = opt_d.clp;
end
end
if opt0.v <= 6 %% convert version 6 to 7
if isfield(opt_d, 'intlinprog')
opt0.intlinprog = opt_d.intlinprog;
end
end
if opt0.v <= 7 %% convert version 7 to 8
opt0.mips.linsolver = opt_d.mips.linsolver;
end
if opt0.v <= 8 %% convert version 8 to 9
opt0.exp.sys_wide_zip_loads = opt_d.exp.sys_wide_zip_loads;
end
opt0.v = v;
end
opt = opt0;
else %% convert from old-style options vector
if DEBUG, fprintf('OPT0 is a old-style options vector\n'); end
opt = mpoption_v2s(opt0);
end
else %% use default options struct as base
if DEBUG, fprintf('no OPT0, starting with default options struct\n'); end
opt = mpoption_default();
end
%% do we have OVERRIDES or NAME/VALUE pairs
ov = [];
if nargin - k == 0 %% looking at last arg, must be OVERRIDES
if isstruct(varargin{k}) %% OVERRIDES provided as struct
if DEBUG, fprintf('OVERRIDES struct\n'); end
ov = varargin{k};
elseif ischar(varargin{k}) %% OVERRIDES provided as file/function name
if DEBUG, fprintf('OVERRIDES file/function name\n'); end
try
ov = feval(varargin{k});
catch
error('mpoption: Unable to load MATPOWER options from ''%s''', varargin{k});
end
if ~isstruct(ov)
error('mpoption: calling ''%s'' did not return a struct', varargin{k});
end
elseif isempty(varargin{k})
return_old_style = 1;
else
error('mpoption: OVERRIDES must be a struct or the name of a function that returns a struct');
end
elseif nargin - k > 0 && mod(nargin-k, 2) %% even number of remaining args
if DEBUG, fprintf('NAME/VALUE pairs override defaults\n'); end
%% process NAME/VALUE pairs
if (have_opt0 && isnumeric(opt0)) ... %% modifying an old-style options vector
|| strcmp(varargin{k}, upper(varargin{k}))
%% this code implies that top-level option fields
%% cannot be all uppercase
if have_opt0
have_old_style_ov = 1;
%% convert pairs to struct
while k < nargin
name = varargin{k};
val = varargin{k+1};
k = k + 2;
ov.(name) = val;
end
else
opt_v = mpoption_old(varargin{:}); %% create modified vector ...
opt = mpoption_v2s(opt_v); %% ... then convert
end
else %% modifying options struct
%% convert pairs to struct
while k < nargin
name = varargin{k};
val = varargin{k+1};
k = k + 2;
c = regexp(name, '([^\.]*)', 'tokens');
s = struct();
for i = 1:length(c)
s(i).type = '.';
s(i).subs = c{i}{1};
end
ov = subsasgn(ov, s, val);
end
end
elseif nargin == 0 || nargin == 1
if DEBUG, fprintf('no OVERRIDES, return default options struct or converted OPT0 vector\n'); end
else
error('mpoption: invalid calling syntax, see ''help mpoption'' to double-check the valid options');
end
%% apply overrides
if ~isempty(ov)
if have_old_style_ov
opt = apply_old_mpoption_overrides(opt, ov);
else
persistent nsc_opt; %% cache this to speed things up
if ~isstruct(nsc_opt)
vf = nested_struct_copy(mpoption_default(), mpoption_info_mips('V'));
vf = nested_struct_copy(vf, mpoption_optional_fields());
ex = struct(...
'name', {...
'cpf.user_callback_args' ...
}, ...
'check', {...
0 ...
}, ...
'copy_mode', {...
'' ...
} ...
);
%% add exceptions for optional packages
opt_pkgs = mpoption_optional_pkgs();
n = length(ex);
for k = 1:length(opt_pkgs)
fname = ['mpoption_info_' opt_pkgs{k}];
if exist(fname, 'file') == 2
opt_ex = feval(fname, 'E');
nex = length(opt_ex);
if ~isempty(opt_ex)
for j = 1:nex
ex(n+j).name = opt_ex(j).name;
end
if isfield(opt_ex, 'check')
for j = 1:nex
ex(n+j).check = opt_ex(j).check;
end
end
if isfield(opt_ex, 'copy_mode')
for j = 1:nex
ex(n+j).copy_mode = opt_ex(j).copy_mode;
end
end
if isfield(opt_ex, 'valid_fields')
for j = 1:nex
ex(n+j).valid_fields = opt_ex(j).valid_fields;
end
end
n = n + nex;
end
end
end
nsc_opt = struct('check', 1, 'valid_fields', vf, 'exceptions', ex);
end
% if have_fcn('catchme')
% try
% opt = nested_struct_copy(opt, ov, nsc_opt);
% catch me
% str = strrep(me.message, 'field', 'option');
% str = strrep(str, 'nested_struct_copy', 'mpoption');
% error(str);
% end
% else
try
opt = nested_struct_copy(opt, ov, nsc_opt);
catch
me = lasterr;
str = strrep(me, 'field', 'option');
str = strrep(str, 'nested_struct_copy', 'mpoption');
error(str);
end
% end
end
end
if return_old_style
opt = mpoption_s2v(opt);
end
%%-------------------------------------------------------------------
function opt = apply_old_mpoption_overrides(opt0, ov)
%
% OPT0 is assumed to already have all of the fields and sub-fields found
% in the default options struct.
%% initialize output
opt = opt0;
errstr = 'mpoption: %g is not a valid value for the old-style ''%s'' option';
fields = fieldnames(ov);
for f = 1:length(fields)
ff = fields{f};
switch ff
case 'PF_ALG'
switch ov.(ff)
case 1
opt.pf.alg = 'NR'; %% Newton's method
case 2
opt.pf.alg = 'FDXB'; %% fast-decoupled (XB version)
case 3
opt.pf.alg = 'FDBX'; %% fast-decoupled (BX version)
case 4
opt.pf.alg = 'GS'; %% Gauss-Seidel
otherwise
error(errstr, ov.(ff), ff);
end
case 'PF_TOL'
opt.pf.tol = ov.(ff);
case 'PF_MAX_IT'
opt.pf.nr.max_it = ov.(ff);
case 'PF_MAX_IT_FD'
opt.pf.fd.max_it = ov.(ff);
case 'PF_MAX_IT_GS'
opt.pf.gs.max_it = ov.(ff);
case 'ENFORCE_Q_LIMS'
opt.pf.enforce_q_lims = ov.(ff);
case 'PF_DC'
switch ov.(ff)
case 0
opt.model = 'AC';
case 1
opt.model = 'DC';
otherwise
error(errstr, ov.(ff), ff);
end
case 'OPF_ALG'
switch ov.(ff)
case 0
opt.opf.ac.solver = 'DEFAULT';
case 500
opt.opf.ac.solver = 'MINOPF';
case 520
opt.opf.ac.solver = 'FMINCON';
case {540, 545}
opt.opf.ac.solver = 'PDIPM';
if ov.(ff) == 545
opt.pdipm.step_control = 1;
else
opt.pdipm.step_control = 0;
end
case 550
opt.opf.ac.solver = 'TRALM';
case {560, 565}
opt.opf.ac.solver = 'MIPS';
if ov.(ff) == 565
opt.mips.step_control = 1;
else
opt.mips.step_control = 0;
end
case 580
opt.opf.ac.solver = 'IPOPT';
case 600
opt.opf.ac.solver = 'KNITRO';
otherwise
error(errstr, ov.(ff), ff);
end
case 'OPF_VIOLATION'
opt.opf.violation = ov.(ff);
case 'CONSTR_TOL_X'
opt.fmincon.tol_x = ov.(ff);
opt.knitro.tol_x = ov.(ff);
case 'CONSTR_TOL_F'
opt.fmincon.tol_f = ov.(ff);
opt.knitro.tol_f = ov.(ff);
case 'CONSTR_MAX_IT'
opt.fmincon.max_it = ov.(ff);
case 'OPF_FLOW_LIM'
switch ov.(ff)
case 0
opt.opf.flow_lim = 'S'; %% apparent power (MVA)
case 1
opt.opf.flow_lim = 'P'; %% real power (MW)
case 2
opt.opf.flow_lim = 'I'; %% current magnitude (MVA @ 1 p.u. voltage)
otherwise
error(errstr, ov.(ff), ff);
end
case 'OPF_IGNORE_ANG_LIM'
opt.opf.ignore_angle_lim = ov.(ff);
case 'OPF_ALG_DC'
switch ov.(ff)
case 0
opt.opf.dc.solver = 'DEFAULT';
case 100
opt.opf.dc.solver = 'BPMPD';
case {200, 250}
opt.opf.dc.solver = 'MIPS';
if ov.(ff) == 250
opt.mips.step_control = 1;
else
opt.mips.step_control = 0;
end
case 300
opt.opf.dc.solver = 'OT'; %% QUADPROG, LINPROG
case 400
opt.opf.dc.solver = 'IPOPT';
case 500
opt.opf.dc.solver = 'CPLEX';
case 600
opt.opf.dc.solver = 'MOSEK';
case 700
opt.opf.dc.solver = 'GUROBI';
otherwise
error(errstr, ov.(ff), ff);
end
case 'VERBOSE'
opt.verbose = ov.(ff);
case 'OUT_ALL'
opt.out.all = ov.(ff);
case 'OUT_SYS_SUM'
opt.out.sys_sum = ov.(ff);
case 'OUT_AREA_SUM'
opt.out.area_sum = ov.(ff);
case 'OUT_BUS'
opt.out.bus = ov.(ff);
case 'OUT_BRANCH'
opt.out.branch = ov.(ff);
case 'OUT_GEN'
opt.out.gen = ov.(ff);
case 'OUT_ALL_LIM'
opt.out.lim.all = ov.(ff);
case 'OUT_V_LIM'
opt.out.lim.v = ov.(ff);
case 'OUT_LINE_LIM'
opt.out.lim.line = ov.(ff);
case 'OUT_PG_LIM'
opt.out.lim.pg = ov.(ff);
case 'OUT_QG_LIM'
opt.out.lim.qg = ov.(ff);
case 'OUT_FORCE'
opt.out.force = ov.(ff);
case 'RETURN_RAW_DER'
opt.opf.return_raw_der = ov.(ff);
case 'FMC_ALG'
opt.fmincon.alg = ov.(ff);
case 'KNITRO_OPT'
opt.knitro.opt = ov.(ff);
case 'IPOPT_OPT'
opt.ipopt.opt = ov.(ff);
case 'MNS_FEASTOL'
opt.minopf.feastol = ov.(ff);
case 'MNS_ROWTOL'
opt.minopf.rowtol = ov.(ff);
case 'MNS_XTOL'
opt.minopf.xtol = ov.(ff);
case 'MNS_MAJDAMP'
opt.minopf.majdamp = ov.(ff);
case 'MNS_MINDAMP'
opt.minopf.mindamp = ov.(ff);
case 'MNS_PENALTY_PARM'
opt.minopf.penalty = ov.(ff);
case 'MNS_MAJOR_IT'
opt.minopf.major_it = ov.(ff);
case 'MNS_MINOR_IT'
opt.minopf.minor_it = ov.(ff);
case 'MNS_MAX_IT'
opt.minopf.max_it = ov.(ff);
case 'MNS_VERBOSITY'
opt.minopf.verbosity = ov.(ff);
case 'MNS_CORE'
opt.minopf.core = ov.(ff);
case 'MNS_SUPBASIC_LIM'
opt.minopf.supbasic_lim = ov.(ff);
case 'MNS_MULT_PRICE'
opt.minopf.mult_price = ov.(ff);
case 'FORCE_PC_EQ_P0'
opt.sopf.force_Pc_eq_P0 = ov.(ff);
case 'PDIPM_FEASTOL'
opt.mips.feastol = ov.(ff);
opt.pdipm.feastol = ov.(ff);
case 'PDIPM_GRADTOL'
opt.mips.gradtol = ov.(ff);
opt.pdipm.gradtol = ov.(ff);
case 'PDIPM_COMPTOL'
opt.mips.comptol = ov.(ff);
opt.pdipm.comptol = ov.(ff);
case 'PDIPM_COSTTOL'
opt.mips.costtol = ov.(ff);
opt.pdipm.costtol = ov.(ff);
case 'PDIPM_MAX_IT'
opt.mips.max_it = ov.(ff);
opt.pdipm.max_it = ov.(ff);
case 'SCPDIPM_RED_IT'
opt.mips.sc.red_it = ov.(ff);
opt.pdipm.sc.red_it = ov.(ff);
case 'TRALM_FEASTOL'
opt.tralm.feastol = ov.(ff);
case 'TRALM_PRIMETOL'
opt.tralm.primaltol = ov.(ff);
case 'TRALM_DUALTOL'
opt.tralm.dualtol = ov.(ff);
case 'TRALM_COSTTOL'
opt.tralm.costtol = ov.(ff);
case 'TRALM_MAJOR_IT'
opt.tralm.major_it = ov.(ff);
case 'TRALM_MINOR_IT'
opt.tralm.minor_it = ov.(ff);
case 'SMOOTHING_RATIO'
opt.pdipm.sc.smooth_ratio = ov.(ff);
opt.tralm.smooth_ratio = ov.(ff);
case 'CPLEX_LPMETHOD'
opt.cplex.lpmethod = ov.(ff);
case 'CPLEX_QPMETHOD'
opt.cplex.qpmethod = ov.(ff);
case 'CPLEX_OPT'
opt.cplex.opt = ov.(ff);
case 'MOSEK_LP_ALG'
opt.mosek.lp_alg = ov.(ff);
case 'MOSEK_MAX_IT'
opt.mosek.max_it = ov.(ff);
case 'MOSEK_GAP_TOL'
opt.mosek.gap_tol = ov.(ff);
case 'MOSEK_MAX_TIME'
opt.mosek.max_time = ov.(ff);
case 'MOSEK_NUM_THREADS'
opt.mosek.num_threads = ov.(ff);
case 'MOSEK_OPT'
opt.mosek.opt = ov.(ff);
case 'GRB_METHOD'
opt.gurobi.method = ov.(ff);
case 'GRB_TIMELIMIT'
opt.gurobi.timelimit = ov.(ff);
case 'GRB_THREADS'
opt.gurobi.threads = ov.(ff);
case 'GRB_OPT'
opt.gurobi.opt = ov.(ff);
otherwise
error('mpoption: ''%s'' is not a valid old-style option name', ff);
end
end
% ov
%%-------------------------------------------------------------------
function opt_s = mpoption_v2s(opt_v)
if DEBUG, fprintf('mpoption_v2s()\n'); end
opt_s = mpoption_default();
errstr = 'mpoption: %g is not a valid value for the old-style ''%s'' option';
switch opt_v(1) %% PF_ALG
case 1
opt_s.pf.alg = 'NR'; %% Newton's method
case 2
opt_s.pf.alg = 'FDXB'; %% fast-decoupled (XB version)
case 3
opt_s.pf.alg = 'FDBX'; %% fast-decoupled (BX version)
case 4
opt_s.pf.alg = 'GS'; %% Gauss-Seidel
otherwise
error(errstr, opt_v(1), 'PF_ALG');
end
opt_s.pf.tol = opt_v(2); %% PF_TOL
opt_s.pf.nr.max_it = opt_v(3); %% PF_MAX_IT
opt_s.pf.fd.max_it = opt_v(4); %% PF_MAX_IT_FD
opt_s.pf.gs.max_it = opt_v(5); %% PF_MAX_IT_GS
opt_s.pf.enforce_q_lims = opt_v(6); %% ENFORCE_Q_LIMS
switch opt_v(10) %% PF_DC
case 0
opt_s.model = 'AC';
case 1
opt_s.model = 'DC';
otherwise
error(errstr, opt_v(10), 'PF_DC');
end
switch opt_v(11) %% OPF_ALG
case 0
opt_s.opf.ac.solver = 'DEFAULT';
case 500
opt_s.opf.ac.solver = 'MINOPF';
case 520
opt_s.opf.ac.solver = 'FMINCON';
case {540, 545}
opt_s.opf.ac.solver = 'PDIPM';
case 550
opt_s.opf.ac.solver = 'TRALM';
case {560, 565}
opt_s.opf.ac.solver = 'MIPS';
case 580
opt_s.opf.ac.solver = 'IPOPT';
case 600
opt_s.opf.ac.solver = 'KNITRO';
otherwise
error(errstr, opt_v(11), 'OPF_ALG');
end
opt_s.opf.violation = opt_v(16); %% OPF_VIOLATION
opt_s.fmincon.tol_x = opt_v(17); %% CONSTR_TOL_X
opt_s.fmincon.tol_f = opt_v(18); %% CONSTR_TOL_F
opt_s.fmincon.max_it = opt_v(19); %% CONSTR_MAX_IT
opt_s.knitro.tol_x = opt_v(17); %% CONSTR_TOL_X
opt_s.knitro.tol_f = opt_v(18); %% CONSTR_TOL_F
switch opt_v(24) %% OPF_FLOW_LIM
case 0
opt_s.opf.flow_lim = 'S'; %% apparent power (MVA)
case 1
opt_s.opf.flow_lim = 'P'; %% real power (MW)
case 2
opt_s.opf.flow_lim = 'I'; %% current magnitude (MVA @ 1 p.u. voltage)
otherwise
error(errstr, opt_v(10), 'PF_DC');
end
opt_s.opf.ignore_angle_lim = opt_v(25); %% OPF_IGNORE_ANG_LIM
switch opt_v(26) %% OPF_ALG_DC
case 0
opt_s.opf.dc.solver = 'DEFAULT';
case 100
opt_s.opf.dc.solver = 'BPMPD';
case {200, 250}
opt_s.opf.dc.solver = 'MIPS';
case 300
opt_s.opf.dc.solver = 'OT'; %% QUADPROG, LINPROG
case 400
opt_s.opf.dc.solver = 'IPOPT';
case 500
opt_s.opf.dc.solver = 'CPLEX';
case 600
opt_s.opf.dc.solver = 'MOSEK';
case 700
opt_s.opf.dc.solver = 'GUROBI';
otherwise
error(errstr, opt_v(26), 'OPF_ALG_DC');
end
opt_s.verbose = opt_v(31); %% VERBOSE
opt_s.out.all = opt_v(32); %% OUT_ALL
opt_s.out.sys_sum = opt_v(33); %% OUT_SYS_SUM
opt_s.out.area_sum = opt_v(34); %% OUT_AREA_SUM
opt_s.out.bus = opt_v(35); %% OUT_BUS
opt_s.out.branch = opt_v(36); %% OUT_BRANCH
opt_s.out.gen = opt_v(37); %% OUT_GEN
opt_s.out.lim.all = opt_v(38); %% OUT_ALL_LIM
opt_s.out.lim.v = opt_v(39); %% OUT_V_LIM
opt_s.out.lim.line = opt_v(40); %% OUT_LINE_LIM
opt_s.out.lim.pg = opt_v(41); %% OUT_PG_LIM
opt_s.out.lim.qg = opt_v(42); %% OUT_QG_LIM
opt_s.out.force = opt_v(44); %% OUT_FORCE
opt_s.opf.return_raw_der = opt_v(52); %% RETURN_RAW_DER
opt_s.fmincon.alg = opt_v(55); %% FMC_ALG
opt_s.knitro.opt = opt_v(58); %% KNITRO_OPT
opt_s.ipopt.opt = opt_v(60); %% IPOPT_OPT
opt_s.minopf.feastol = opt_v(61); %% MNS_FEASTOL
opt_s.minopf.rowtol = opt_v(62); %% MNS_ROWTOL
opt_s.minopf.xtol = opt_v(63); %% MNS_XTOL
opt_s.minopf.majdamp = opt_v(64); %% MNS_MAJDAMP
opt_s.minopf.mindamp = opt_v(65); %% MNS_MINDAMP
opt_s.minopf.penalty = opt_v(66); %% MNS_PENALTY_PARM
opt_s.minopf.major_it = opt_v(67); %% MNS_MAJOR_IT
opt_s.minopf.minor_it = opt_v(68); %% MNS_MINOR_IT
opt_s.minopf.max_it = opt_v(69); %% MNS_MAX_IT
opt_s.minopf.verbosity = opt_v(70); %% MNS_VERBOSITY
opt_s.minopf.core = opt_v(71); %% MNS_CORE
opt_s.minopf.supbasic_lim = opt_v(72); %% MNS_SUPBASIC_LIM
opt_s.minopf.mult_price = opt_v(73); %% MNS_MULT_PRICE
opt_s.sopf.force_Pc_eq_P0 = opt_v(80); %% FORCE_PC_EQ_P0, for c3sopf
if (opt_v(11) == 565 && opt_v(10) == 0) || (opt_v(26) == 250 && opt_v(10) == 1)
opt_s.mips.step_control = 1;
end
opt_s.mips.feastol = opt_v(81); %% PDIPM_FEASTOL
opt_s.mips.gradtol = opt_v(82); %% PDIPM_GRADTOL
opt_s.mips.comptol = opt_v(83); %% PDIPM_COMPTOL
opt_s.mips.costtol = opt_v(84); %% PDIPM_COSTTOL
opt_s.mips.max_it = opt_v(85); %% PDIPM_MAX_IT
opt_s.mips.sc.red_it = opt_v(86); %% SCPDIPM_RED_IT
opt_s.pdipm.feastol = opt_v(81); %% PDIPM_FEASTOL
opt_s.pdipm.gradtol = opt_v(82); %% PDIPM_GRADTOL
opt_s.pdipm.comptol = opt_v(83); %% PDIPM_COMPTOL
opt_s.pdipm.costtol = opt_v(84); %% PDIPM_COSTTOL
opt_s.pdipm.max_it = opt_v(85); %% PDIPM_MAX_IT
opt_s.pdipm.sc.red_it = opt_v(86); %% SCPDIPM_RED_IT
opt_s.pdipm.sc.smooth_ratio = opt_v(93); %% SMOOTHING_RATIO
if opt_v(11) == 545 && opt_v(10) == 0
opt_s.pdipm.step_control = 1;
end
opt_s.tralm.feastol = opt_v(87); %% TRALM_FEASTOL
opt_s.tralm.primaltol = opt_v(88); %% TRALM_PRIMETOL
opt_s.tralm.dualtol = opt_v(89); %% TRALM_DUALTOL
opt_s.tralm.costtol = opt_v(90); %% TRALM_COSTTOL
opt_s.tralm.major_it = opt_v(91); %% TRALM_MAJOR_IT
opt_s.tralm.minor_it = opt_v(92); %% TRALM_MINOR_IT
opt_s.tralm.smooth_ratio = opt_v(93); %% SMOOTHING_RATIO
opt_s.cplex.lpmethod = opt_v(95); %% CPLEX_LPMETHOD
opt_s.cplex.qpmethod = opt_v(96); %% CPLEX_QPMETHOD
opt_s.cplex.opt = opt_v(97); %% CPLEX_OPT
opt_s.mosek.lp_alg = opt_v(111); %% MOSEK_LP_ALG
opt_s.mosek.max_it = opt_v(112); %% MOSEK_MAX_IT
opt_s.mosek.gap_tol = opt_v(113); %% MOSEK_GAP_TOL
opt_s.mosek.max_time = opt_v(114); %% MOSEK_MAX_TIME
opt_s.mosek.num_threads = opt_v(115); %% MOSEK_NUM_THREADS
opt_s.mosek.opt = opt_v(116); %% MOSEK_OPT
opt_s.gurobi.method = opt_v(121); %% GRB_METHOD
opt_s.gurobi.timelimit = opt_v(122); %% GRB_TIMELIMIT
opt_s.gurobi.threads = opt_v(123); %% GRB_THREADS
opt_s.gurobi.opt = opt_v(124); %% GRB_OPT
function opt_v = mpoption_s2v(opt_s)
if DEBUG, fprintf('mpoption_s2v()\n'); end
%% PF_ALG
old = mpoption_old;
switch upper(opt_s.pf.alg)
case 'NR'
PF_ALG = 1;
case 'FDXB'
PF_ALG = 2;
case 'FDBX'
PF_ALG = 3;
case 'GS'
PF_ALG = 4;
end
%% PF_DC
if strcmp(upper(opt_s.model), 'DC')
PF_DC = 1;
else
PF_DC = 0;
end
%% OPF_ALG
switch upper(opt_s.opf.ac.solver)
case 'DEFAULT'
OPF_ALG = 0;
case 'MINOPF'
OPF_ALG = 500;
case 'FMINCON'
OPF_ALG = 520;
case 'PDIPM'
if isfield(opt_s, 'pdipm') && opt_s.pdipm.step_control
OPF_ALG = 545;
else
OPF_ALG = 540;
end
case 'TRALM'
OPF_ALG = 550;
case 'MIPS'
if opt_s.mips.step_control
OPF_ALG = 565;
else
OPF_ALG = 560;
end
case 'IPOPT'
OPF_ALG = 580;
case 'KNITRO'
OPF_ALG = 600;
end
%% FMINCON, Knitro tol_x, tol_f, max_it
if strcmp(upper(opt_s.opf.ac.solver), 'KNITRO') && isfield(opt_s, 'knitro')
CONSTR_TOL_X = opt_s.knitro.tol_x;
CONSTR_TOL_F = opt_s.knitro.tol_f;
elseif isfield(opt_s, 'fmincon')
CONSTR_TOL_X = opt_s.fmincon.tol_x;
CONSTR_TOL_F = opt_s.fmincon.tol_f;
else
CONSTR_TOL_X = old(17);
CONSTR_TOL_F = old(18);
end
if isfield(opt_s, 'fmincon')
CONSTR_MAX_IT = opt_s.fmincon.max_it;
FMC_ALG = opt_s.fmincon.alg;
else
CONSTR_MAX_IT = old(19);
FMC_ALG = old(55);
end
%% OPF_FLOW_LIM
switch upper(opt_s.opf.flow_lim)
case 'S'
OPF_FLOW_LIM = 0;
case 'P'
OPF_FLOW_LIM = 1;
case 'I'
OPF_FLOW_LIM = 2;
end
%% OPF_ALG_DC
switch upper(opt_s.opf.dc.solver)
case 'DEFAULT'
OPF_ALG_DC = 0;
case 'BPMPD'
OPF_ALG_DC = 100;
case 'MIPS'
if opt_s.mips.step_control
OPF_ALG_DC = 250;
else
OPF_ALG_DC = 200;
end
case 'OT'
OPF_ALG_DC = 300;
case 'IPOPT'
OPF_ALG_DC = 400;
case 'CPLEX'
OPF_ALG_DC = 500;
case 'MOSEK'
OPF_ALG_DC = 600;
case 'GUROBI'
OPF_ALG_DC = 700;
end
%% KNITRO_OPT
if isfield(opt_s, 'knitro')
KNITRO_OPT = opt_s.knitro.opt;
else
KNITRO_OPT = old(58);
end
%% IPOPT_OPT
if isfield(opt_s, 'ipopt')
IPOPT_OPT = opt_s.ipopt.opt;
else
IPOPT_OPT = old(58);
end
%% MINOPF options
if isfield(opt_s, 'minopf')
MINOPF_OPTS = [
opt_s.minopf.feastol; %% 61 - MNS_FEASTOL
opt_s.minopf.rowtol; %% 62 - MNS_ROWTOL
opt_s.minopf.xtol; %% 63 - MNS_XTOL
opt_s.minopf.majdamp; %% 64 - MNS_MAJDAMP
opt_s.minopf.mindamp; %% 65 - MNS_MINDAMP
opt_s.minopf.penalty; %% 66 - MNS_PENALTY_PARM
opt_s.minopf.major_it; %% 67 - MNS_MAJOR_IT
opt_s.minopf.minor_it; %% 68 - MNS_MINOR_IT
opt_s.minopf.max_it; %% 69 - MNS_MAX_IT
opt_s.minopf.verbosity; %% 70 - MNS_VERBOSITY
opt_s.minopf.core; %% 71 - MNS_CORE
opt_s.minopf.supbasic_lim; %% 72 - MNS_SUPBASIC_LIM
opt_s.minopf.mult_price;%% 73 - MNS_MULT_PRICE
];
else
MINOPF_OPTS = old(61:73);
end
%% FORCE_PC_EQ_P0
if isfield(opt_s, 'sopf') && isfield(opt_s.sopf, 'force_Pc_eq_P0')
FORCE_PC_EQ_P0 = opt_s.sopf.force_Pc_eq_P0;
else
FORCE_PC_EQ_P0 = 0;
end
%% PDIPM options
if isfield(opt_s, 'pdipm')
PDIPM_OPTS = [
opt_s.pdipm.feastol; %% 81 - PDIPM_FEASTOL
opt_s.pdipm.gradtol; %% 82 - PDIPM_GRADTOL
opt_s.pdipm.comptol; %% 83 - PDIPM_COMPTOL
opt_s.pdipm.costtol; %% 84 - PDIPM_COSTTOL
opt_s.pdipm.max_it; %% 85 - PDIPM_MAX_IT
opt_s.pdipm.sc.red_it; %% 86 - SCPDIPM_RED_IT
];
else
PDIPM_OPTS = old(81:86);
end
%% TRALM options
if isfield(opt_s, 'tralm')
TRALM_OPTS = [
opt_s.tralm.feastol; %% 87 - TRALM_FEASTOL
opt_s.tralm.primaltol; %% 88 - TRALM_PRIMETOL
opt_s.tralm.dualtol; %% 89 - TRALM_DUALTOL
opt_s.tralm.costtol; %% 90 - TRALM_COSTTOL
opt_s.tralm.major_it; %% 91 - TRALM_MAJOR_IT
opt_s.tralm.minor_it; %% 92 - TRALM_MINOR_IT
];
else
TRALM_OPTS = old(87:92);
end
%% SMOOTHING_RATIO
if strcmp(upper(opt_s.opf.ac.solver), 'TRALM') && isfield(opt_s, 'tralm')
SMOOTHING_RATIO = opt_s.tralm.smooth_ratio;
elseif isfield(opt_s, 'pdipm')
SMOOTHING_RATIO = opt_s.pdipm.sc.smooth_ratio;
else
SMOOTHING_RATIO = old(93);
end
%% CPLEX options
if isfield(opt_s, 'cplex')
CPLEX_OPTS = [
opt_s.cplex.lpmethod; %% 95 - CPLEX_LPMETHOD
opt_s.cplex.qpmethod; %% 96 - CPLEX_QPMETHOD
opt_s.cplex.opt; %% 97 - CPLEX_OPT
];
else
CPLEX_OPTS = old(95:97);
end
%% MOSEK options
if isfield(opt_s, 'mosek')
MOSEK_OPTS = [
opt_s.mosek.lp_alg; %% 111 - MOSEK_LP_ALG
opt_s.mosek.max_it; %% 112 - MOSEK_MAX_IT
opt_s.mosek.gap_tol; %% 113 - MOSEK_GAP_TOL
opt_s.mosek.max_time; %% 114 - MOSEK_MAX_TIME
opt_s.mosek.num_threads;%% 115 - MOSEK_NUM_THREADS
opt_s.mosek.opt; %% 116 - MOSEK_OPT
];
else
MOSEK_OPTS = old(111:116);
end
%% Gurobi options
if isfield(opt_s, 'gurobi')
GUROBI_OPTS = [
opt_s.gurobi.method; %% 121 - GRB_METHOD
opt_s.gurobi.timelimit; %% 122 - GRB_TIMELIMIT
opt_s.gurobi.threads; %% 123 - GRB_THREADS
opt_s.gurobi.opt; %% 124 - GRB_OPT
];
else
GUROBI_OPTS = old(121:124);
end
opt_v = [
%% power flow options
PF_ALG; %% 1 - PF_ALG
opt_s.pf.tol; %% 2 - PF_TOL
opt_s.pf.nr.max_it; %% 3 - PF_MAX_IT
opt_s.pf.fd.max_it; %% 4 - PF_MAX_IT_FD
opt_s.pf.gs.max_it; %% 5 - PF_MAX_IT_GS
opt_s.pf.enforce_q_lims;%% 6 - ENFORCE_Q_LIMS
0; %% 7 - RESERVED7
0; %% 8 - RESERVED8
0; %% 9 - RESERVED9
PF_DC; %% 10 - PF_DC
%% OPF options
OPF_ALG; %% 11 - OPF_ALG
0; %% 12 - RESERVED12 (was OPF_ALG_POLY = 100)
0; %% 13 - RESERVED13 (was OPF_ALG_PWL = 200)
0; %% 14 - RESERVED14 (was OPF_POLY2PWL_PTS = 10)
0; %% 15 - OPF_NEQ (removed)
opt_s.opf.violation; %% 16 - OPF_VIOLATION
CONSTR_TOL_X; %% 17 - CONSTR_TOL_X
CONSTR_TOL_F; %% 18 - CONSTR_TOL_F
CONSTR_MAX_IT; %% 19 - CONSTR_MAX_IT
old(20); %% 20 - LPC_TOL_GRAD (removed)
old(21); %% 21 - LPC_TOL_X (removed)
old(22); %% 22 - LPC_MAX_IT (removed)
old(23); %% 23 - LPC_MAX_RESTART (removed)
OPF_FLOW_LIM; %% 24 - OPF_FLOW_LIM
opt_s.opf.ignore_angle_lim; %% 25 - OPF_IGNORE_ANG_LIM
OPF_ALG_DC; %% 26 - OPF_ALG_DC
0; %% 27 - RESERVED27
0; %% 28 - RESERVED28
0; %% 29 - RESERVED29
0; %% 30 - RESERVED30
%% output options
opt_s.verbose; %% 31 - VERBOSE
opt_s.out.all; %% 32 - OUT_ALL
opt_s.out.sys_sum; %% 33 - OUT_SYS_SUM
opt_s.out.area_sum; %% 34 - OUT_AREA_SUM
opt_s.out.bus; %% 35 - OUT_BUS
opt_s.out.branch; %% 36 - OUT_BRANCH
opt_s.out.gen; %% 37 - OUT_GEN
opt_s.out.lim.all; %% 38 - OUT_ALL_LIM
opt_s.out.lim.v; %% 39 - OUT_V_LIM
opt_s.out.lim.line; %% 40 - OUT_LINE_LIM
opt_s.out.lim.pg; %% 41 - OUT_PG_LIM
opt_s.out.lim.qg; %% 42 - OUT_QG_LIM
0; %% 43 - RESERVED43 (was OUT_RAW)
opt_s.out.force; %% 44 - OUT_FORCE
0; %% 45 - RESERVED45
0; %% 46 - RESERVED46
0; %% 47 - RESERVED47
0; %% 48 - RESERVED48
0; %% 49 - RESERVED49
0; %% 50 - RESERVED50
%% other options
old(51); %% 51 - SPARSE_QP (removed)
opt_s.opf.return_raw_der; %% 52 - RETURN_RAW_DER
0; %% 53 - RESERVED53
0; %% 54 - RESERVED54
FMC_ALG; %% 55 - FMC_ALG
0; %% 56 - RESERVED56
0; %% 57 - RESERVED57
KNITRO_OPT; %% 58 - KNITRO_OPT
0; %% 59 - RESERVED59
IPOPT_OPT; %% 60 - IPOPT_OPT
%% MINOPF options
MINOPF_OPTS; %% 61-73 - MNS_FEASTOL-MNS_MULT_PRICE
0; %% 74 - RESERVED74
0; %% 75 - RESERVED75
0; %% 76 - RESERVED76
0; %% 77 - RESERVED77
0; %% 78 - RESERVED78
0; %% 79 - RESERVED79
FORCE_PC_EQ_P0; %% 80 - FORCE_PC_EQ_P0, for c3sopf
%% MIPS, PDIPM, SC-PDIPM, and TRALM options
PDIPM_OPTS; %% 81-86 - PDIPM_FEASTOL-SCPDIPM_RED_IT
TRALM_OPTS; %% 87-92 - TRALM_FEASTOL-TRALM_MINOR_IT
SMOOTHING_RATIO; %% 93 - SMOOTHING_RATIO
0; %% 94 - RESERVED94
%% CPLEX options
CPLEX_OPTS; %% 95-97 - CPLEX_LPMETHOD-CPLEX_OPT
0; %% 98 - RESERVED98
0; %% 99 - RESERVED99
0; %% 100 - RESERVED100
0; %% 101 - RESERVED101
0; %% 102 - RESERVED102
0; %% 103 - RESERVED103
0; %% 104 - RESERVED104
0; %% 105 - RESERVED105
0; %% 106 - RESERVED106
0; %% 107 - RESERVED107
0; %% 108 - RESERVED108
0; %% 109 - RESERVED109
0; %% 110 - RESERVED110
%% MOSEK options
MOSEK_OPTS; %% 111-116 - MOSEK_LP_ALG-MOSEK_OPT
0; %% 117 - RESERVED117
0; %% 118 - RESERVED118
0; %% 119 - RESERVED119
0; %% 120 - RESERVED120
%% Gurobi options
GUROBI_OPTS; %% 121-124 - GRB_METHOD-GRB_OPT
];
%%-------------------------------------------------------------------
function optt = mpoption_default()
if DEBUG, fprintf('mpoption_default()\n'); end
persistent opt; %% cache this for speed
if ~isstruct(opt)
opt = struct(...
'v', mpoption_version, ... %% version
'model', 'AC', ...
'pf', struct(...
'alg', 'NR', ...
'tol', 1e-8, ...
'nr', struct(...
'max_it', 10 ), ...
'fd', struct(...
'max_it', 30 ), ...
'gs', struct(...
'max_it', 1000 ), ...
'enforce_q_lims', 0 ), ...
'cpf', struct(...
'parameterization', 3, ...
'stop_at', 'NOSE', ... %% 'NOSE', <lam val>, 'FULL'
'step', 0.05, ...
'adapt_step', 0, ...
'error_tol', 1e-3, ...
'step_min', 1e-4, ...
'step_max', 0.2, ...
'plot', struct(...
'level', 0, ...
'bus', [] ), ...
'user_callback', '', ...
'user_callback_args', struct() ), ...
'opf', struct(...
'ac', struct(...
'solver', 'DEFAULT' ), ...
'dc', struct(...
'solver', 'DEFAULT' ), ...
'violation', 5e-6, ...
'flow_lim', 'S', ...
'ignore_angle_lim', 0, ...
'init_from_mpc', -1, ...
'return_raw_der', 0 ), ...
'verbose', 1, ...
'out', struct(...
'all', -1, ...
'sys_sum', 1, ...
'area_sum', 0, ...
'bus', 1, ...
'branch', 1, ...
'gen', 0, ...
'lim', struct(...
'all', -1, ...
'v', 1, ...
'line', 1, ...
'pg', 1, ...
'qg', 1 ), ...
'force', 0, ...
'suppress_detail', -1 ), ...
'mips', struct(... %% see mpoption_info_mips() for optional fields
'step_control', 0, ...
'linsolver', '', ...
'feastol', 0, ...
'gradtol', 1e-6, ...
'comptol', 1e-6, ...
'costtol', 1e-6, ...
'max_it', 150, ...
'sc', struct(...
'red_it', 20 )), ...
'exp', struct(... %% experimental options
'sys_wide_zip_loads', struct(...
'pw', [], ...
'qw', [] )) ...
);
opt_pkgs = mpoption_optional_pkgs();
for k = 1:length(opt_pkgs)
fname = ['mpoption_info_' opt_pkgs{k}];
if exist(fname, 'file') == 2
opt = nested_struct_copy(opt, feval(fname, 'D'));
end
end
end
optt = opt;
%%-------------------------------------------------------------------
function optt = mpoption_optional_fields()
if DEBUG, fprintf('mpoption_optional_fields()\n'); end
persistent opt; %% cache this for speed
if ~isstruct(opt)
opt_pkgs = mpoption_optional_pkgs();
opt = struct;
for k = 1:length(opt_pkgs)
fname = ['mpoption_info_' opt_pkgs{k}];
if exist(fname, 'file') == 2
opt = nested_struct_copy(opt, feval(fname, 'V'));
end
end
end
optt = opt;
%% globals
%%-------------------------------------------------------------------
function v = mpoption_version
v = 10; %% version number of MATPOWER options struct
%% (must be incremented every time structure is updated)
%% v1 - first version based on struct (MATPOWER 5.0b1)
%% v2 - added 'linprog' and 'quadprog' fields
%% v3 - (forgot to increment v) added 'out.suppress_detail'
%% field
%% v4 - (forgot to increment v) MIPS 1.1, added optional
%% fields to 'mips' options: xi, sigma, z0, alpha_min,
%% rho_min, rho_max, mu_threshold and max_stepsize
%% v5 - (forgot to increment v) added 'opf.init_from_mpc'
%% field (MATPOWER 5.0)
%% v6 - added 'clp' field
%% v7 - added 'intlinprog' field
%% v8 - MIPS 1.2, added 'linsolver' field to
%% 'mips' options
%% v9 - added 'exp' for experimental fields, specifically
%% 'sys_wide_zip_loads.pw', 'sys_wide_zip_loads.qw'
%% v10 - added 'most' field
%%-------------------------------------------------------------------
function db_level = DEBUG
db_level = 0;
%%-------------------------------------------------------------------
function pkgs = mpoption_optional_pkgs()
pkgs = {...
'clp', 'cplex', 'fmincon', 'gurobi', 'glpk', 'intlinprog', 'ipopt', ...
'knitro', 'linprog', 'minopf', 'most', 'mosek', 'quadprog', 'sdp_pf', ...
'sopf', 'tspopf', 'yalmip' ...
};
function [bus, gen, branch] = pfsoln(baseMVA, bus0, gen0, branch0, Ybus, Yf, Yt, V, ref, pv, pq, mpopt)
%PFSOLN Updates bus, gen, branch data structures to match power flow soln.
% [BUS, GEN, BRANCH] = PFSOLN(BASEMVA, BUS0, GEN0, BRANCH0, ...
% YBUS, YF, YT, V, REF, PV, PQ, MPOPT)
% MATPOWER
% Copyright (c) 1996-2016 by Power System Engineering Research Center (PSERC)
% by Ray Zimmerman, PSERC Cornell
%
% This file is part of MATPOWER.
% Covered by the 3-clause BSD License (see LICENSE file for details).
% See http://www.pserc.cornell.edu/matpower/ for more info.
%% define named indices into bus, gen, branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
%% default options
if nargin < 12
mpopt = mpoption();
end
%% initialize return values
bus = bus0;
gen = gen0;
branch = branch0;
%%----- update bus voltages -----
bus(:, VM) = abs(V);
bus(:, VA) = angle(V) * 180 / pi;
%%----- update Qg for gens at PV/slack buses and Pg for slack bus(es) -----
%% generator info
on = find(gen(:, GEN_STATUS) > 0 & ... %% which generators are on?
bus(gen(:, GEN_BUS), BUS_TYPE) ~= PQ); %% ... and not at PQ buses
off = find(gen(:, GEN_STATUS) <= 0); %% which generators are off?
gbus = gen(on, GEN_BUS); %% what buses are they at?
%% compute total injected bus powers
Sbus = V(gbus) .* conj(Ybus(gbus, :) * V);
%% update Qg for generators at PV/slack buses
gen(off, QG) = zeros(length(off), 1); %% zero out off-line Qg
%% don't touch the ones at PQ buses
[Pd_gbus, Qd_gbus] = total_load(bus(gbus, :), [], 'bus', [], mpopt);
gen(on, QG) = imag(Sbus) * baseMVA + Qd_gbus; %% inj Q + local Qd
%% ... at this point any buses with more than one generator will have
%% the total Q dispatch for the bus assigned to each generator. This
%% must be split between them. We do it first equally, then in proportion
%% to the reactive range of the generator.
if length(on) > 1
%% build connection matrix, element i, j is 1 if gen on(i) at bus j is ON
nb = size(bus, 1);
ngon = size(on, 1);
Cg = sparse((1:ngon)', gbus, ones(ngon, 1), ngon, nb);
%% divide Qg by number of generators at the bus to distribute equally
ngg = Cg * sum(Cg)'; %% ngon x 1, number of gens at this gen's bus
gen(on, QG) = gen(on, QG) ./ ngg;
%% divide proportionally
Cmin = sparse((1:ngon)', gbus, gen(on, QMIN), ngon, nb);
Cmax = sparse((1:ngon)', gbus, gen(on, QMAX), ngon, nb);
Qg_tot = Cg' * gen(on, QG); %% nb x 1 vector of total Qg at each bus
Qg_min = sum(Cmin)'; %% nb x 1 vector of min total Qg at each bus
Qg_max = sum(Cmax)'; %% nb x 1 vector of max total Qg at each bus
ig = find(Cg * Qg_min == Cg * Qg_max); %% gens at buses with Qg range = 0
Qg_save = gen(on(ig), QG);
gen(on, QG) = gen(on, QMIN) + ...
(Cg * ((Qg_tot - Qg_min)./(Qg_max - Qg_min + eps))) .* ...
(gen(on, QMAX) - gen(on, QMIN)); %% ^ avoid div by 0
gen(on(ig), QG) = Qg_save;
end %% (terms are mult by 0 anyway)
%% update Pg for slack gen(s)
for k = 1:length(ref)
refgen = find(gbus == ref(k)); %% which is(are) the reference gen(s)?
Pd_refk = total_load(bus(ref(k), :), [], 'bus', [], mpopt);
gen(on(refgen(1)), PG) = real(Sbus(refgen(1))) * baseMVA + Pd_refk; %% inj P + local Pd
if length(refgen) > 1 %% more than one generator at this ref bus
%% subtract off what is generated by other gens at this bus
gen(on(refgen(1)), PG) = gen(on(refgen(1)), PG) ...
- sum(gen(on(refgen(2:length(refgen))), PG));
end
end
%%----- update/compute branch power flows -----
out = find(branch(:, BR_STATUS) == 0); %% out-of-service branches
br = find(branch(:, BR_STATUS)); %% in-service branches
Sf = V(branch(br, F_BUS)) .* conj(Yf(br, :) * V) * baseMVA; %% complex power at "from" bus
St = V(branch(br, T_BUS)) .* conj(Yt(br, :) * V) * baseMVA; %% complex power injected at "to" bus
branch(br, [PF, QF, PT, QT]) = [real(Sf) imag(Sf) real(St) imag(St)];
branch(out, [PF, QF, PT, QT]) = zeros(length(out), 4);
function [Pd, Qd] = total_load(bus, gen, load_zone, opt, mpopt)
%TOTAL_LOAD Returns vector of total load in each load zone.
% PD = TOTAL_LOAD(MPC)
% PD = TOTAL_LOAD(MPC, LOAD_ZONE)
% PD = TOTAL_LOAD(MPC, LOAD_ZONE, OPT)
% PD = TOTAL_LOAD(MPC, LOAD_ZONE, OPT, MPOPT)
% PD = TOTAL_LOAD(BUS)
% PD = TOTAL_LOAD(BUS, GEN)
% PD = TOTAL_LOAD(BUS, GEN, LOAD_ZONE)
% PD = TOTAL_LOAD(BUS, GEN, LOAD_ZONE, OPT)
% PD = TOTAL_LOAD(BUS, GEN, LOAD_ZONE, OPT, MPOPT)
% [PD, QD] = TOTAL_LOAD(...) returns both active and reative power
% demand for each zone.
%
% MPC - standard MATPOWER case struct
%
% BUS - standard BUS matrix with nb rows, where the fixed active
% and reactive loads are specified in columns PD and QD
%
% GEN - (optional) standard GEN matrix with ng rows, where the
% dispatchable loads are specified by columns PG, QG, PMIN,
% QMIN and QMAX (in rows for which ISLOAD(GEN) returns true).
% If GEN is empty, it assumes there are no dispatchable loads.
%
% LOAD_ZONE - (optional) nb element vector where the value of
% each element is either zero or the index of the load zone
% to which the corresponding bus belongs. If LOAD_ZONE(b) = k
% then the loads at bus b will added to the values of PD(k) and
% QD(k). If LOAD_ZONE is empty, the default is defined as the areas
% specified in the BUS matrix, i.e. LOAD_ZONE = BUS(:, BUS_AREA)
% and load will have dimension = MAX(BUS(:, BUS_AREA)). LOAD_ZONE
% can also take the following string values:
% 'all' - use a single zone for the entire system (return scalar)
% 'area' - use LOAD_ZONE = BUS(:, BUS_AREA), same as default
% 'bus' - use a different zone for each bus (i.e. to compute
% final values of bus-wise loads, including voltage dependent
% fixed loads and or dispatchable loads)
%
% OPT - (optional) option struct, with the following fields:
% 'type' - string specifying types of loads to include, default
% is 'BOTH' if GEN is provided, otherwise 'FIXED'
% 'FIXED' : sum only fixed loads
% 'DISPATCHABLE' : sum only dispatchable loads
% 'BOTH' : sum both fixed and dispatchable loads
% 'nominal' - 1 : use nominal load for dispatchable loads
% 0 : (default) use actual realized load for
% dispatchable loads
%
% For backward compatibility with MATPOWER 4.x, OPT can also
% take the form of a string, with the same options as OPT.type above.
% In this case, again for backward compatibility, it is the "nominal"
% load that is computed for dispatchable loads, not the actual
% realized load. Using a string for OPT is deprecated and
% will be removed in a future version.
%
% MPOPT - (optional) MATPOWER options struct, which may specify
% a voltage dependent (ZIP) load model for fixed loads
%
% Examples:
% Return the total active load for each area as defined in BUS_AREA.
%
% Pd = total_load(bus);
%
% Return total active and reactive load, fixed and dispatchable, for
% entire system.
%
% [Pd, Qd] = total_load(bus, gen, 'all');
%
% Return the total of the nominal dispatchable loads at buses 10-20.
%
% load_zone = zeros(nb, 1);
% load_zone(10:20) = 1;
% opt = struct('type', 'DISPATCHABLE', 'nominal', 1);
% Pd = total_load(mpc, load_zone, opt)
%
% See also SCALE_LOAD.
% MATPOWER
% Copyright (c) 2004-2016 by Power System Engineering Research Center (PSERC)
% by Ray Zimmerman, PSERC Cornell
%
% This file is part of MATPOWER.
% Covered by the 3-clause BSD License (see LICENSE file for details).
% See http://www.pserc.cornell.edu/matpower/ for more info.
%% define constants
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
%% purposely being backward compatible with older MATPOWER
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, ...
PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen;
%%----- process inputs -----
if nargin < 5
mpopt = [];
if nargin < 4
opt = [];
if nargin < 3
load_zone = [];
if nargin < 2
gen = [];
end
end
end
end
if isstruct(bus) %% shift input arguments
mpopt = opt;
opt = load_zone;
load_zone = gen;
mpc = bus;
bus = mpc.bus;
gen = mpc.gen;
end
nb = size(bus, 1); %% number of buses
%% default options
if ischar(opt) %% convert old WHICH_TYPE string option to struct
opt = struct('type', opt, 'nominal', 1);
else
if ~isfield(opt, 'type')
if isempty(gen)
opt.type = 'FIXED';
else
opt.type = 'BOTH';
end
end
if ~isfield(opt, 'nominal')
opt.nominal = 0;
end
end
switch upper(opt.type(1))
case {'F', 'D', 'B'}
%% OK
otherwise
error('total_load: OPT.type should be ''FIXED'', ''DISPATCHABLE'' or ''BOTH''');
end
want_Q = (nargout > 1);
want_fixed = (opt.type(1) == 'B' || opt.type(1) == 'F');
want_disp = (opt.type(1) == 'B' || opt.type(1) == 'D');
%% initialize load_zone
if ischar(load_zone)
if strcmp(lower(load_zone), 'bus')
load_zone = (1:nb)'; %% each bus is its own zone
elseif strcmp(lower(load_zone), 'all')
load_zone = ones(nb, 1); %% make a single zone of all buses
elseif strcmp(lower(load_zone), 'area')
load_zone = bus(:, BUS_AREA); %% use areas defined in bus data as zones
end
elseif isempty(load_zone)
load_zone = bus(:, BUS_AREA); %% use areas defined in bus data as zones
end
nz = max(load_zone); %% number of load zones
%% fixed load at each bus, & initialize dispatchable
if want_fixed
Sd = makeSdzip(1, bus, mpopt);
Vm = bus(:, VM);
Sbusd = Sd.p + Sd.i .* Vm + Sd.z .* Vm.^2;
Pdf = real(Sbusd); %% real power
if want_Q
Qdf = imag(Sbusd); %% reactive power
end
else
Pdf = zeros(nb, 1); %% real power
if want_Q
Qdf = zeros(nb, 1); %% reactive power
end
end
%% dispatchable load at each bus
if want_disp %% need dispatchable
ng = size(gen, 1);
is_ld = isload(gen) & gen(:, GEN_STATUS) > 0;
ld = find(is_ld);
%% create map of external bus numbers to bus indices
i2e = bus(:, BUS_I);
e2i = sparse(max(i2e), 1);
e2i(i2e) = (1:nb)';
Cld = sparse(e2i(gen(:, GEN_BUS)), (1:ng)', is_ld, nb, ng);
if opt.nominal %% use nominal power
Pdd = -Cld * gen(:, PMIN); %% real power
if want_Q
Q = zeros(ng, 1);
Q(ld) = (gen(ld, QMIN) == 0) .* gen(ld, QMAX) + ...
(gen(ld, QMAX) == 0) .* gen(ld, QMIN);
Qdd = -Cld * Q; %% reactive power
end
else %% use realized actual power dispatch
Pdd = -Cld * gen(:, PG); %% real power
if want_Q
Qdd = -Cld * gen(:, QG); %% reactive power
end
end
else
Pdd = zeros(nb, 1);
if want_Q
Qdd = zeros(nb, 1);
end
end
%% compute load sums
if nz == nb && all(load_zone == (1:nb)'); %% individual buses
Pd = (Pdf + Pdd) .* (bus(:, BUS_TYPE) ~= NONE);
if want_Q
Qd = (Qdf + Qdd) .* (bus(:, BUS_TYPE) ~= NONE);
end
else
Pd = zeros(nz, 1);
if want_Q
Qd = zeros(nz, 1);
end
for k = 1:nz
idx = find( load_zone == k & bus(:, BUS_TYPE) ~= NONE);
Pd(k) = sum(Pdf(idx)) + sum(Pdd(idx));
if want_Q
Qd(k) = sum(Qdf(idx)) + sum(Qdd(idx));
end
end
end
142
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
