-
Notifications
You must be signed in to change notification settings - Fork 149
/
opf.m
305 lines (290 loc) · 12.6 KB
/
opf.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
function [busout, genout, branchout, f, success, info, et, g, jac, xr, pimul] = ...
opf(varargin)
% opf - Solves an optimal power flow.
% ::
%
% [RESULTS, SUCCESS] = OPF(MPC, MPOPT)
%
% Returns either a RESULTS struct and an optional SUCCESS flag, or individual
% data matrices, the objective function value and a SUCCESS flag. In the
% latter case, there are additional optional return values. See Examples
% below for the possible calling syntax options.
%
% Examples:
% Output argument options:
%
% results = opf(...)
% [results, success] = opf(...)
% [bus, gen, branch, f, success] = opf(...)
% [bus, gen, branch, f, success, info, et, g, jac, xr, pimul] = opf(...)
%
% Input arguments options:
%
% opf(mpc)
% opf(mpc, mpopt)
% opf(mpc, userfcn, mpopt)
% opf(mpc, A, l, u)
% opf(mpc, A, l, u, mpopt)
% opf(mpc, A, l, u, mpopt, N, fparm, H, Cw)
% opf(mpc, A, l, u, mpopt, N, fparm, H, Cw, z0, zl, zu)
%
% opf(baseMVA, bus, gen, branch, areas, gencost)
% opf(baseMVA, bus, gen, branch, areas, gencost, mpopt)
% opf(baseMVA, bus, gen, branch, areas, gencost, userfcn, mpopt)
% opf(baseMVA, bus, gen, branch, areas, gencost, A, l, u)
% opf(baseMVA, bus, gen, branch, areas, gencost, A, l, u, mpopt)
% opf(baseMVA, bus, gen, branch, areas, gencost, A, l, u, ...
% mpopt, N, fparm, H, Cw)
% opf(baseMVA, bus, gen, branch, areas, gencost, A, l, u, ...
% mpopt, N, fparm, H, Cw, z0, zl, zu)
%
% The data for the problem can be specified in one of three ways:
% (1) a string (mpc) containing the file name of a MATPOWER case
% which defines the data matrices baseMVA, bus, gen, branch, and
% gencost (areas is not used at all, it is only included for
% backward compatibility of the API).
% (2) a struct (mpc) containing the data matrices as fields.
% (3) the individual data matrices themselves.
%
% The optional user parameters for user constraints (A, l, u), user costs
% (N, fparm, H, Cw), user variable initializer (z0), and user variable
% limits (zl, zu) can also be specified as fields in a case struct,
% either passed in directly or defined in a case file referenced by name.
%
% When specified, A, l, u represent additional linear constraints on the
% optimization variables, l <= A*[x; z] <= u. If the user specifies an A
% matrix that has more columns than the number of "x" (OPF) variables,
% then there are extra linearly constrained "z" variables. For an
% explanation of the formulation used and instructions for forming the
% A matrix, see the manual.
%
% A generalized cost on all variables can be applied if input arguments
% N, fparm, H and Cw are specified. First, a linear transformation
% of the optimization variables is defined by means of r = N * [x; z].
% Then, to each element of r a function is applied as encoded in the
% fparm matrix (see manual). If the resulting vector is named w,
% then H and Cw define a quadratic cost on w: (1/2)*w'*H*w + Cw * w .
% H and N should be sparse matrices and H should also be symmetric.
%
% The optional mpopt vector specifies MATPOWER options. If the OPF
% algorithm is not explicitly set in the options MATPOWER will use
% the default solver, based on a primal-dual interior point method.
% For the AC OPF this is opf.ac.solver = 'MIPS', unless the TSPOPF optional
% package is installed, in which case the default is 'PDIPM'. For the
% DC OPF, the default is opf.dc.solver = 'MIPS'. See MPOPTION for
% more details on the available OPF solvers and other OPF options
% and their default values.
%
% The solved case is returned either in a single results struct (described
% below) or in the individual data matrices, bus, gen and branch. Also
% returned are the final objective function value (f) and a flag which is
% true if the algorithm was successful in finding a solution (success).
% Additional optional return values are an algorithm specific return status
% (info), elapsed time in seconds (et), the constraint vector (g), the
% Jacobian matrix (jac), and the vector of variables (xr) as well
% as the constraint multipliers (pimul).
%
% The single results struct is a MATPOWER case struct (mpc) with the
% usual baseMVA, bus, branch, gen, gencost fields, along with the
% following additional fields:
%
% .order see 'help ext2int' for details of this field
% .et elapsed time in seconds for solving OPF
% .success 1 if solver converged successfully, 0 otherwise
% .om OPF model object, see 'help opf_model'
% .x final value of optimization variables (internal order)
% .f final objective function value
% .mu shadow prices on ...
% .var
% .l lower bounds on variables
% .u upper bounds on variables
% .nln
% .l lower bounds on nonlinear constraints
% .u upper bounds on nonlinear constraints
% .lin
% .l lower bounds on linear constraints
% .u upper bounds on linear constraints
% .raw raw solver output in form returned by MINOS, and more
% .xr final value of optimization variables
% .pimul constraint multipliers
% .info solver specific termination code
% .output solver specific output information
% .alg algorithm code of solver used
% .g (optional) constraint values
% .dg (optional) constraint 1st derivatives
% .df (optional) obj fun 1st derivatives (not yet implemented)
% .d2f (optional) obj fun 2nd derivatives (not yet implemented)
% .var
% .val optimization variable values, by named block
% .Va voltage angles
% .Vm voltage magnitudes (AC only)
% .Pg real power injections
% .Qg reactive power injections (AC only)
% .y constrained cost variable (only if have pwl costs)
% (other) any user defined variable blocks
% .mu variable bound shadow prices, by named block
% .l lower bound shadow prices
% .Va, Vm, Pg, Qg, y, (other)
% .u upper bound shadow prices
% .Va, Vm, Pg, Qg, y, (other)
% .nle (AC only)
% .lambda shadow prices on nonlinear equality constraints,
% by named block
% .Pmis real power mismatch equations
% .Qmis reactive power mismatch equations
% (other) use defined constraints
% .nli (AC only)
% .mu shadow prices on nonlinear inequality constraints,
% by named block
% .Sf flow limits at "from" end of branches
% .St flow limits at "to" end of branches
% (other) use defined constraints
% .lin
% .mu shadow prices on linear constraints, by named block
% .l lower bounds
% .Pmis real power mistmatch equations (DC only)
% .Pf flow limits at "from" end of branches (DC only)
% .Pt flow limits at "to" end of branches (DC only)
% .PQh upper portion of gen PQ-capability curve (AC only)
% .PQl lower portion of gen PQ-capability curve (AC only)
% .vl constant power factor constraint for loads (AC only)
% .ycon basin constraints for CCV for pwl costs
% (other) any user defined constraint blocks
% .u upper bounds
% .Pmis, Pf, Pt, PQh, PQl, vl, ycon, (other)
% .cost user defined cost values, by named block
%
% See also runopf, dcopf, uopf, caseformat.
% MATPOWER
% Copyright (c) 1996-2024, Power Systems Engineering Research Center (PSERC)
% by Ray Zimmerman, PSERC Cornell
% and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
%
% This file is part of MATPOWER.
% Covered by the 3-clause BSD License (see LICENSE file for details).
% See https://matpower.org for more info.
%%----- initialization -----
t0 = tic; %% start timer
%% define named indices into data 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;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
%% process input arguments
[mpc, mpopt] = opf_args(varargin{:});
%% if 'opf.ac.solver' not set, choose MIPS
if strcmp(upper(mpopt.opf.ac.solver), 'DEFAULT')
mpopt = mpoption(mpopt, 'opf.ac.solver', 'MIPS'); %% MIPS
end
%% initialize state with power flow solution, if requested
if mpopt.opf.start == 3
mpopt_pf = mpoption(mpopt, 'out.all', 0, 'verbose', max(0, mpopt.verbose-1));
if mpopt.verbose
fprintf('Running power flow to initialize OPF.\n');
end
rpf = runpf(mpc, mpopt_pf);
if rpf.success
mpc = rpf; %% or should I just copy Va, Vm, Pg, Qg?
end
end
%% add zero columns to bus, gen, branch for multipliers, etc if needed
nb = size(mpc.bus, 1); %% number of buses
nl = size(mpc.branch, 1); %% number of branches
ng = size(mpc.gen, 1); %% number of dispatchable injections
if size(mpc.bus,2) < MU_VMIN
mpc.bus = [mpc.bus zeros(nb, MU_VMIN-size(mpc.bus,2)) ];
end
if size(mpc.gen,2) < MU_QMIN
mpc.gen = [ mpc.gen zeros(ng, MU_QMIN-size(mpc.gen,2)) ];
end
if size(mpc.branch,2) < MU_ANGMAX
mpc.branch = [ mpc.branch zeros(nl, MU_ANGMAX-size(mpc.branch,2)) ];
end
%%----- convert to internal numbering, remove out-of-service stuff -----
mpc = ext2int(mpc, mpopt);
%%----- construct OPF model object -----
%% use MP-Core?
have_mp_core = have_feature('mp_core');
use_mp_core = 0;
if have_mp_core && ~mpopt.exp.use_legacy_core
dc = strcmp(upper(mpopt.model), 'DC');
alg = upper(mpopt.opf.ac.solver);
if dc || ~(strcmp(alg, 'MINOPF') || strcmp(alg, 'PDIPM') || ...
strcmp(alg, 'TRALM') || strcmp(alg, 'SDPOPF'))
use_mp_core = 1;
end
end
if ~use_mp_core
om = opf_setup(mpc, mpopt);
end
%%----- execute the OPF -----
if nargout > 7
mpopt.opf.return_raw_der = 1;
end
if use_mp_core
task_class = @mp.task_opf_legacy; %% set default task class
%% get and apply extensions
if isfield(mpopt.exp, 'mpx') && ~isempty(mpopt.exp.mpx)
mpx = mpopt.exp.mpx;
if ~iscell(mpx)
mpx = { mpx };
end
else
mpx = {};
end
for k = 1:length(mpx)
task_class = mpx{k}.task_class(task_class, mpopt);
end
%% create and run task
opf = task_class();
opf.run(mpc, mpopt, mpx);
[results, success, raw] = opf.legacy_post_run(mpopt);
else
if ~isempty(mpc.bus)
[results, success, raw] = opf_execute(om, mpopt);
else
results = mpc;
success = 0;
raw.output.message = 'MATPOWER case contains no connected buses';
if mpopt.verbose
fprintf('OPF not valid : %s\n', raw.output.message);
end
end
results.success = success; %% make success available to subsequent callbacks
end
%%----- revert to original ordering, including out-of-service stuff -----
results = int2ext(results);
%% 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 MU_PMAX MU_PMIN]) = 0;
end
if ~isempty(results.order.branch.status.off)
results.branch(results.order.branch.status.off, [PF QF PT QT MU_SF MU_ST MU_ANGMIN MU_ANGMAX]) = 0;
end
%%----- finish preparing output -----
et = toc(t0); %% compute elapsed time
if nargout > 0
if nargout <= 2
results.et = et;
results.raw = raw;
busout = results;
genout = success;
else
[busout, genout, branchout, f, info, xr, pimul] = deal(results.bus, ...
results.gen, results.branch, results.f, raw.info, raw.xr, raw.pimul);
if isfield(results, 'g')
g = results.g;
end
if isfield(results, 'dg')
jac = results.dg;
end
end
elseif success
results.et = et;
printpf(results, 1, mpopt);
end