0001 function [res, suc] = ...
0002 runcpf(basecasedata, targetcasedata, mpopt, fname, solvedcase)
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0088 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0089 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0090 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0091 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0092 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0093 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0094 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0095
0096
0097
0098 if nargin < 5
0099 solvedcase = '';
0100 if nargin < 4
0101 fname = '';
0102 if nargin < 3
0103 mpopt = mpoption;
0104 if nargin < 2
0105 targetcasedata = 'case9target';
0106 if nargin < 1
0107 basecasedata = 'case9';
0108 end
0109 end
0110 end
0111 end
0112 end
0113
0114
0115 step = mpopt.cpf.step;
0116 parameterization = mpopt.cpf.parameterization;
0117 adapt_step = mpopt.cpf.adapt_step;
0118 cb_args = mpopt.cpf.user_callback_args;
0119
0120
0121 callback_names = {'cpf_default_callback'};
0122 if ~isempty(mpopt.cpf.user_callback)
0123 if iscell(mpopt.cpf.user_callback)
0124 callback_names = {callback_names{:}, mpopt.cpf.user_callback{:}};
0125 else
0126 callback_names = {callback_names{:}, mpopt.cpf.user_callback};
0127 end
0128 end
0129 callbacks = cellfun(@str2func, callback_names, 'UniformOutput', false);
0130
0131
0132 mpcbase = loadcase(basecasedata);
0133 nb = size(mpcbase.bus, 1);
0134
0135
0136 if size(mpcbase.branch,2) < QT
0137 mpcbase.branch = [ mpcbase.branch zeros(size(mpcbase.branch, 1), QT-size(mpcbase.branch,2)) ];
0138 end
0139
0140
0141 mpcbase = ext2int(mpcbase);
0142 [baseMVAb, busb, genb, branchb] = deal(mpcbase.baseMVA, mpcbase.bus, mpcbase.gen, mpcbase.branch);
0143
0144
0145 [ref, pv, pq] = bustypes(busb, genb);
0146
0147
0148 onb = find(genb(:, GEN_STATUS) > 0);
0149 gbusb = genb(onb, GEN_BUS);
0150
0151
0152 mpctarget = loadcase(targetcasedata);
0153
0154
0155 if size(mpctarget.branch,2) < QT
0156 mpctarget.branch = [ mpctarget.branch zeros(size(mpctarget.branch, 1), QT-size(mpctarget.branch,2)) ];
0157 end
0158
0159
0160 mpctarget = ext2int(mpctarget);
0161 [baseMVAt, bust, gent, brancht] = deal(mpctarget.baseMVA, mpctarget.bus, mpctarget.gen, mpctarget.branch);
0162
0163
0164
0165
0166
0167 ont = find(gent(:, GEN_STATUS) > 0);
0168 gbust = gent(ont, GEN_BUS);
0169
0170
0171 t0 = clock;
0172 if mpopt.verbose > 0
0173 v = mpver('all');
0174 fprintf('\nMATPOWER Version %s, %s', v.Version, v.Date);
0175 fprintf(' -- AC Continuation Power Flow\n');
0176 end
0177
0178
0179 V0 = busb(:, VM) .* exp(sqrt(-1) * pi/180 * busb(:, VA));
0180 vcb = ones(size(V0));
0181 vcb(pq) = 0;
0182 k = find(vcb(gbusb));
0183 V0(gbusb(k)) = genb(onb(k), VG) ./ abs(V0(gbusb(k))).* V0(gbusb(k));
0184
0185
0186 [Ybus, Yf, Yt] = makeYbus(baseMVAb, busb, branchb);
0187
0188
0189 Sbusb = makeSbus(baseMVAb, busb, genb);
0190
0191 Sbust = makeSbus(baseMVAt, bust, gent);
0192
0193
0194 Sxfr = Sbust - Sbusb;
0195
0196
0197 if mpopt.verbose > 2
0198 mpopt_pf = mpoption(mpopt, 'verbose', max(0, mpopt.verbose-1));
0199 else
0200 mpopt_pf = mpoption(mpopt, 'verbose', max(0, mpopt.verbose-2));
0201 end
0202 lam = 0;
0203 [V, success, iterations] = newtonpf(Ybus, Sbusb, V0, ref, pv, pq, mpopt_pf);
0204 if mpopt.verbose > 2
0205 fprintf('step %3d : lambda = %6.3f\n', 0, 0);
0206 elseif mpopt.verbose > 1
0207 fprintf('step %3d : lambda = %6.3f, %2d Newton steps\n', 0, 0, iterations);
0208 end
0209
0210 lamprv = lam;
0211 Vprv = V;
0212 continuation = 1;
0213 cont_steps = 0;
0214
0215
0216 cb_data = struct( ...
0217 'mpc_base', mpcbase, ...
0218 'mpc_target', mpctarget, ...
0219 'Sxfr', Sxfr, ...
0220 'Ybus', Ybus, ...
0221 'Yf', Yf, ...
0222 'Yt', Yt, ...
0223 'ref', ref, ...
0224 'pv', pv, ...
0225 'pq', pq, ...
0226 'mpopt', mpopt );
0227 cb_state = struct();
0228
0229
0230 for k = 1:length(callbacks)
0231 cb_state = callbacks{k}(cont_steps, V, lam, V, lam, ...
0232 cb_data, cb_state, cb_args);
0233 end
0234
0235 if norm(Sxfr) == 0
0236 if mpopt.verbose
0237 fprintf('base case and target case have identical load and generation\n');
0238 end
0239 continuation = 0;
0240 V0 = V; lam0 = lam;
0241 end
0242
0243
0244 z = zeros(2*length(V)+1,1);
0245 z(end,1) = 1.0;
0246 while(continuation)
0247 cont_steps = cont_steps + 1;
0248
0249 [V0, lam0, z] = cpf_predictor(V, lam, Ybus, Sxfr, pv, pq, step, z, ...
0250 Vprv, lamprv, parameterization);
0251
0252
0253 Vprv = V;
0254 lamprv = lam;
0255
0256
0257 [V, success, i, lam] = cpf_corrector(Ybus, Sbusb, V0, ref, pv, pq, ...
0258 lam0, Sxfr, Vprv, lamprv, z, step, parameterization, mpopt_pf);
0259 if ~success
0260 continuation = 0;
0261 if mpopt.verbose
0262 fprintf('step %3d : lambda = %6.3f, corrector did not converge in %d iterations\n', cont_steps, lam, i);
0263 end
0264 break;
0265 end
0266 if mpopt.verbose > 2
0267 fprintf('step %3d : lambda = %6.3f\n', cont_steps, lam);
0268 elseif mpopt.verbose > 1
0269 fprintf('step %3d : lambda = %6.3f, %2d corrector Newton steps\n', cont_steps, lam, i);
0270 end
0271
0272
0273 for k = 1:length(callbacks)
0274 cb_state = callbacks{k}(cont_steps, V, lam, V0, lam0, ...
0275 cb_data, cb_state, cb_args);
0276 end
0277
0278 if ischar(mpopt.cpf.stop_at)
0279 if strcmp(upper(mpopt.cpf.stop_at), 'FULL')
0280 if abs(lam) < 1e-8
0281
0282 if mpopt.verbose
0283 fprintf('\nTraced full continuation curve in %d continuation steps\n',cont_steps);
0284 end
0285 continuation = 0;
0286 elseif lam < lamprv && lam - step < 0
0287 step = lam;
0288 parameterization = 1;
0289 adapt_step = 0;
0290 end
0291 else
0292 if lam < lamprv
0293 if mpopt.verbose
0294 fprintf('\nReached steady state loading limit in %d continuation steps\n',cont_steps);
0295 end
0296 continuation = 0;
0297 end
0298 end
0299 else
0300 if lam < lamprv
0301 if mpopt.verbose
0302 fprintf('\nReached steady state loading limit in %d continuation steps\n', cont_steps);
0303 end
0304 continuation = 0;
0305 elseif abs(mpopt.cpf.stop_at - lam) < 1e-8
0306 if mpopt.verbose
0307 fprintf('\nReached desired lambda %3.2f in %d continuation steps\n', ...
0308 mpopt.cpf.stop_at, cont_steps);
0309 end
0310 continuation = 0;
0311 elseif lam + step > mpopt.cpf.stop_at
0312 step = mpopt.cpf.stop_at - lam;
0313 parameterization = 1;
0314 adapt_step = 0;
0315 end
0316 end
0317
0318 if adapt_step && continuation
0319
0320 cpf_error = norm([angle(V(pq));abs(V([pv;pq]));lam] - [angle(V0(pq));abs(V0([pv;pq]));lam0],inf);
0321 if cpf_error < mpopt.cpf.error_tol
0322
0323 step = step*mpopt.cpf.error_tol/cpf_error;
0324 if step > mpopt.cpf.step_max
0325 step = mpopt.cpf.step_max;
0326 end
0327 else
0328
0329 step = step*mpopt.cpf.error_tol/cpf_error;
0330 if step < mpopt.cpf.step_min
0331 step = mpopt.cpf.step_min;
0332 end
0333 end
0334 end
0335 end
0336
0337
0338 if success
0339 cpf_results = struct();
0340 for k = 1:length(callbacks)
0341 [cb_state, cpf_results] = callbacks{k}(cont_steps, V, lam, V0, lam0, ...
0342 cb_data, cb_state, cb_args, cpf_results);
0343 end
0344 else
0345 cpf_results.iterations = i;
0346 end
0347
0348
0349
0350 bust(:,PD) = busb(:,PD) + lam*(bust(:,PD) - busb(:,PD));
0351 bust(:,QD) = busb(:,QD) + lam*(bust(:,QD) - busb(:,QD));
0352 gent(:,PG) = genb(:,PG) + lam*(gent(:,PG) - genb(:,PG));
0353
0354
0355 [bust, gent, brancht] = pfsoln(baseMVAt, bust, gent, brancht, Ybus, Yf, Yt, V, ref, pv, pq);
0356
0357 mpctarget.et = etime(clock, t0);
0358 mpctarget.success = success;
0359
0360
0361
0362 [mpctarget.bus, mpctarget.gen, mpctarget.branch] = deal(bust, gent, brancht);
0363 if success
0364 n = cpf_results.iterations + 1;
0365 cpf_results.V_p = i2e_data(mpctarget, cpf_results.V_p, NaN(nb,n), 'bus', 1);
0366 cpf_results.V_c = i2e_data(mpctarget, cpf_results.V_c, NaN(nb,n), 'bus', 1);
0367 end
0368 results = int2ext(mpctarget);
0369 results.cpf = cpf_results;
0370
0371
0372 if ~isempty(results.order.gen.status.off)
0373 results.gen(results.order.gen.status.off, [PG QG]) = 0;
0374 end
0375 if ~isempty(results.order.branch.status.off)
0376 results.branch(results.order.branch.status.off, [PF QF PT QT]) = 0;
0377 end
0378
0379 if fname
0380 [fd, msg] = fopen(fname, 'at');
0381 if fd == -1
0382 error(msg);
0383 else
0384 if mpopt.out.all == 0
0385 printpf(results, fd, mpoption(mpopt, 'out.all', -1));
0386 else
0387 printpf(results, fd, mpopt);
0388 end
0389 fclose(fd);
0390 end
0391 end
0392 printpf(results, 1, mpopt);
0393
0394
0395 if solvedcase
0396 savecase(solvedcase, results);
0397 end
0398
0399 if nargout
0400 res = results;
0401 if nargout > 1
0402 suc = success;
0403 end
0404
0405 end