-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathrunpf_complete.m
5536 lines (5181 loc) · 233 KB
/
runpf_complete.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
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
% MATPOWER CODES COMNBINED TOGETHER IN ONE FILE TO AVOID DEPENDENCIES
function [MVAbase, bus, gen, branch, success, et,J,Ybus] = ...
runpf_complete(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_complete;
% results = runpf_complete(casedata);
% results = runpf_complete(casedata, mpopt);
% results = runpf_complete(casedata, mpopt, fname);
% results = runpf_complete(casedata, mpopt, fname, solvedcase);
% [results, success] = runpf_complete(...);
%
% 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, Power Systems 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 = 'case30'; %% default data file is 'case9.m'
end
end
end
end
% mpopt.pf.enforce_q_lims=1;
%% 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;
its = 0; %% total iterations
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);
its = 1;
%% 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,J] = 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
its = its + iterations;
%% 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
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) && ...
(isempty(mx) || isempty(mn))
%% all remaining PV/REF gens are violating AND all are
%% violating same limit (all violating Qmin or all Qmax)
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;
mpc.iterations = its;
%%----- 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 involved in this
function varargout = deal(varargin)
%DEAL Deal inputs to outputs.
% [A,B,C,...] = DEAL(X,Y,Z,...) simply matches up the input and
% output lists. It is the same as A=X, B=Y, C=Z, ...
% [A,B,C,...] = DEAL(X) copies the single input to all
% the requested outputs. It is the same as A=X, B=X, C=X, ...
%
% DEAL is most useful when used with cell arrays and structures
% via comma separated list expansion. Here are some useful
% constructions:
% [S.FIELD] = DEAL(X) sets all the fields with the name FIELD
% in the structure array S to the value X. If S doesn't
% exist, use [S(1:M).FIELD] = DEAL(X);
% [X{:}] = DEAL(A.FIELD) copies the values of the field with
% name FIELD to the cell array X. If X doesn't exist,
% use [X{1:M}] = DEAL(A.FIELD).
% [A,B,C,...] = DEAL(X{:}) copies the contents of the cell
% array X to the separate variables A,B,C,...
% [A,B,C,...] = DEAL(S.FIELD) copies the contents of the fields
% with the name FIELD to separate variables A,B,C,...
%
% Examples:
% sys = {rand(3) ones(3,1) eye(3) zeros(3,1)};
% [a,b,c,d] = deal(sys{:});
%
% direc = dir; filenames = {};
% [filenames{1:length(direc),1}] = deal(direc.name);
%
% See also LISTS, PAREN.
% Copyright 1984-2005 The MathWorks, Inc.
if nargin==1,
varargout = varargin(ones(1,nargout));
else
if nargout ~= nargin
error(message('MATLAB:deal:narginNargoutMismatch'))
end
varargout = varargin;
end
function [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
%IDX_BUS Defines constants for named column indices to bus matrix.
% Example:
%
% [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;
%
% Some examples of usage, after defining the constants using the line above,
% are:
%
% Pd = bus(4, PD); % get the real power demand at bus 4
% bus(:, VMIN) = 0.95; % set the min voltage magnitude to 0.95 at all buses
%
% The index, name and meaning of each column of the bus matrix is given
% below:
%
% columns 1-13 must be included in input matrix (in case file)
% 1 BUS_I bus number (positive integer)
% 2 BUS_TYPE bus type (1 = PQ, 2 = PV, 3 = ref, 4 = isolated)
% 3 PD Pd, real power demand (MW)
% 4 QD Qd, reactive power demand (MVAr)
% 5 GS Gs, shunt conductance (MW demanded at V = 1.0 p.u.)
% 6 BS Bs, shunt susceptance (MVAr injected at V = 1.0 p.u.)
% 7 BUS_AREA area number, (positive integer)
% 8 VM Vm, voltage magnitude (p.u.)
% 9 VA Va, voltage angle (degrees)
% 10 BASE_KV baseKV, base voltage (kV)
% 11 ZONE zone, loss zone (positive integer)
% 12 VMAX maxVm, maximum voltage magnitude (p.u.)
% 13 VMIN minVm, minimum voltage magnitude (p.u.)
%
% columns 14-17 are added to matrix after OPF solution
% they are typically not present in the input matrix
% (assume OPF objective function has units, u)
% 14 LAM_P Lagrange multiplier on real power mismatch (u/MW)
% 15 LAM_Q Lagrange multiplier on reactive power mismatch (u/MVAr)
% 16 MU_VMAX Kuhn-Tucker multiplier on upper voltage limit (u/p.u.)
% 17 MU_VMIN Kuhn-Tucker multiplier on lower voltage limit (u/p.u.)
%
% additional constants, used to assign/compare values in the BUS_TYPE column
% 1 PQ PQ bus
% 2 PV PV bus
% 3 REF reference bus
% 4 NONE isolated bus
%
% See also DEFINE_CONSTANTS.
% MATPOWER
% Copyright (c) 1996-2016, Power Systems 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 bus types
PQ = 1;
PV = 2;
REF = 3;
NONE = 4;
%% define the indices
BUS_I = 1; %% bus number (1 to 29997)
BUS_TYPE = 2; %% bus type (1 - PQ bus, 2 - PV bus, 3 - reference bus, 4 - isolated bus)
PD = 3; %% Pd, real power demand (MW)
QD = 4; %% Qd, reactive power demand (MVAr)
GS = 5; %% Gs, shunt conductance (MW at V = 1.0 p.u.)
BS = 6; %% Bs, shunt susceptance (MVAr at V = 1.0 p.u.)
BUS_AREA = 7; %% area number, 1-100
VM = 8; %% Vm, voltage magnitude (p.u.)
VA = 9; %% Va, voltage angle (degrees)
BASE_KV = 10; %% baseKV, base voltage (kV)
ZONE = 11; %% zone, loss zone (1-999)
VMAX = 12; %% maxVm, maximum voltage magnitude (p.u.) (not in PTI format)
VMIN = 13; %% minVm, minimum voltage magnitude (p.u.) (not in PTI format)
%% included in opf solution, not necessarily in input
%% assume objective function has units, u
LAM_P = 14; %% Lagrange multiplier on real power mismatch (u/MW)
LAM_Q = 15; %% Lagrange multiplier on reactive power mismatch (u/MVAr)
MU_VMAX = 16; %% Kuhn-Tucker multiplier on upper voltage limit (u/p.u.)
MU_VMIN = 17; %% Kuhn-Tucker multiplier on lower voltage limit (u/p.u.)
function [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
%IDX_GEN Defines constants for named column indices to gen matrix.
% Example:
%
% [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;
%
% Some examples of usage, after defining the constants using the line above,
% are:
%
% Pg = gen(4, PG); % get the real power output of generator 4
% gen(:, PMIN) = 0; % set to zero the minimum real power limit of all gens
%
% The index, name and meaning of each column of the gen matrix is given
% below:
%
% columns 1-21 must be included in input matrix (in case file)
% 1 GEN_BUS bus number
% 2 PG Pg, real power output (MW)
% 3 QG Qg, reactive power output (MVAr)
% 4 QMAX Qmax, maximum reactive power output (MVAr)
% 5 QMIN Qmin, minimum reactive power output (MVAr)
% 6 VG Vg, voltage magnitude setpoint (p.u.)
% 7 MBASE mBase, total MVA base of machine, defaults to baseMVA
% 8 GEN_STATUS status, > 0 - in service, <= 0 - out of service
% 9 PMAX Pmax, maximum real power output (MW)
% 10 PMIN Pmin, minimum real power output (MW)
% 11 PC1 Pc1, lower real power output of PQ capability curve (MW)
% 12 PC2 Pc2, upper real power output of PQ capability curve (MW)
% 13 QC1MIN Qc1min, minimum reactive power output at Pc1 (MVAr)
% 14 QC1MAX Qc1max, maximum reactive power output at Pc1 (MVAr)
% 15 QC2MIN Qc2min, minimum reactive power output at Pc2 (MVAr)
% 16 QC2MAX Qc2max, maximum reactive power output at Pc2 (MVAr)
% 17 RAMP_AGC ramp rate for load following/AGC (MW/min)
% 18 RAMP_10 ramp rate for 10 minute reserves (MW)
% 19 RAMP_30 ramp rate for 30 minute reserves (MW)
% 20 RAMP_Q ramp rate for reactive power (2 sec timescale) (MVAr/min)
% 21 APF area participation factor
%
% columns 22-25 are added to matrix after OPF solution
% they are typically not present in the input matrix
% (assume OPF objective function has units, u)
% 22 MU_PMAX Kuhn-Tucker multiplier on upper Pg limit (u/MW)
% 23 MU_PMIN Kuhn-Tucker multiplier on lower Pg limit (u/MW)
% 24 MU_QMAX Kuhn-Tucker multiplier on upper Qg limit (u/MVAr)
% 25 MU_QMIN Kuhn-Tucker multiplier on lower Qg limit (u/MVAr)
%
% See also DEFINE_CONSTANTS.
% MATPOWER
% Copyright (c) 1996-2016, Power Systems 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 the indices
GEN_BUS = 1; %% bus number
PG = 2; %% Pg, real power output (MW)
QG = 3; %% Qg, reactive power output (MVAr)
QMAX = 4; %% Qmax, maximum reactive power output at Pmin (MVAr)
QMIN = 5; %% Qmin, minimum reactive power output at Pmin (MVAr)
VG = 6; %% Vg, voltage magnitude setpoint (p.u.)
MBASE = 7; %% mBase, total MVA base of this machine, defaults to baseMVA
GEN_STATUS = 8; %% status, 1 - machine in service, 0 - machine out of service
PMAX = 9; %% Pmax, maximum real power output (MW)
PMIN = 10; %% Pmin, minimum real power output (MW)
PC1 = 11; %% Pc1, lower real power output of PQ capability curve (MW)
PC2 = 12; %% Pc2, upper real power output of PQ capability curve (MW)
QC1MIN = 13; %% Qc1min, minimum reactive power output at Pc1 (MVAr)
QC1MAX = 14; %% Qc1max, maximum reactive power output at Pc1 (MVAr)
QC2MIN = 15; %% Qc2min, minimum reactive power output at Pc2 (MVAr)
QC2MAX = 16; %% Qc2max, maximum reactive power output at Pc2 (MVAr)
RAMP_AGC = 17; %% ramp rate for load following/AGC (MW/min)
RAMP_10 = 18; %% ramp rate for 10 minute reserves (MW)
RAMP_30 = 19; %% ramp rate for 30 minute reserves (MW)
RAMP_Q = 20; %% ramp rate for reactive power (2 sec timescale) (MVAr/min)
APF = 21; %% area participation factor
%% included in opf solution, not necessarily in input
%% assume objective function has units, u
MU_PMAX = 22; %% Kuhn-Tucker multiplier on upper Pg limit (u/MW)
MU_PMIN = 23; %% Kuhn-Tucker multiplier on lower Pg limit (u/MW)
MU_QMAX = 24; %% Kuhn-Tucker multiplier on upper Qg limit (u/MVAr)
MU_QMIN = 25; %% Kuhn-Tucker multiplier on lower Qg limit (u/MVAr)
%% Note: When a generator's PQ capability curve is not simply a box and the
%% upper Qg limit is binding, the multiplier on this constraint is split into
%% it's P and Q components and combined with the appropriate MU_Pxxx and
%% MU_Qxxx values. Likewise for the lower Q limits.
function [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
%IDX_BRCH Defines constants for named column indices to branch matrix.
% Example:
%
% [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;
%
% Some examples of usage, after defining the constants using the line above,
% are:
%
% branch(4, BR_STATUS) = 0; % take branch 4 out of service
% Ploss = branch(:, PF) + branch(:, PT); % compute real power loss vector
%
% The index, name and meaning of each column of the branch matrix is given
% below:
%
% columns 1-11 must be included in input matrix (in case file)
% 1 F_BUS f, from bus number
% 2 T_BUS t, to bus number
% 3 BR_R r, resistance (p.u.)
% 4 BR_X x, reactance (p.u.)
% 5 BR_B b, total line charging susceptance (p.u.)
% 6 RATE_A rateA, MVA rating A (long term rating)
% 7 RATE_B rateB, MVA rating B (short term rating)
% 8 RATE_C rateC, MVA rating C (emergency rating)
% 9 TAP ratio, transformer off nominal turns ratio
% 10 SHIFT angle, transformer phase shift angle (degrees)
% 11 BR_STATUS initial branch status, 1 - in service, 0 - out of service
% 12 ANGMIN minimum angle difference, angle(Vf) - angle(Vt) (degrees)
% 13 ANGMAX maximum angle difference, angle(Vf) - angle(Vt) (degrees)
% (The voltage angle difference is taken to be unbounded below
% if ANGMIN < -360 and unbounded above if ANGMAX > 360.
% If both parameters are zero, it is unconstrained.)
%
% columns 14-17 are added to matrix after power flow or OPF solution
% they are typically not present in the input matrix
% 14 PF real power injected at "from" bus end (MW)
% 15 QF reactive power injected at "from" bus end (MVAr)
% 16 PT real power injected at "to" bus end (MW)
% 17 QT reactive power injected at "to" bus end (MVAr)
%
% columns 18-21 are added to matrix after OPF solution
% they are typically not present in the input matrix
% (assume OPF objective function has units, u)
% 18 MU_SF Kuhn-Tucker multiplier on MVA limit at "from" bus (u/MVA)
% 19 MU_ST Kuhn-Tucker multiplier on MVA limit at "to" bus (u/MVA)
% 20 MU_ANGMIN Kuhn-Tucker multiplier lower angle difference limit (u/degree)
% 21 MU_ANGMAX Kuhn-Tucker multiplier upper angle difference limit (u/degree)
%
% See also DEFINE_CONSTANTS.
% MATPOWER
% Copyright (c) 1996-2016, Power Systems 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 the indices
F_BUS = 1; %% f, from bus number
T_BUS = 2; %% t, to bus number
BR_R = 3; %% r, resistance (p.u.)
BR_X = 4; %% x, reactance (p.u.)
BR_B = 5; %% b, total line charging susceptance (p.u.)
RATE_A = 6; %% rateA, MVA rating A (long term rating)
RATE_B = 7; %% rateB, MVA rating B (short term rating)
RATE_C = 8; %% rateC, MVA rating C (emergency rating)
TAP = 9; %% ratio, transformer off nominal turns ratio
SHIFT = 10; %% angle, transformer phase shift angle (degrees)
BR_STATUS = 11; %% initial branch status, 1 - in service, 0 - out of service
ANGMIN = 12; %% minimum angle difference, angle(Vf) - angle(Vt) (degrees)
ANGMAX = 13; %% maximum angle difference, angle(Vf) - angle(Vt) (degrees)
%% included in power flow solution, not necessarily in input
PF = 14; %% real power injected at "from" bus end (MW) (not in PTI format)
QF = 15; %% reactive power injected at "from" bus end (MVAr) (not in PTI format)
PT = 16; %% real power injected at "to" bus end (MW) (not in PTI format)
QT = 17; %% reactive power injected at "to" bus end (MVAr) (not in PTI format)
%% included in opf solution, not necessarily in input
%% assume objective function has units, u
MU_SF = 18; %% Kuhn-Tucker multiplier on MVA limit at "from" bus (u/MVA)
MU_ST = 19; %% Kuhn-Tucker multiplier on MVA limit at "to" bus (u/MVA)
MU_ANGMIN = 20; %% Kuhn-Tucker multiplier lower angle difference limit (u/degree)
MU_ANGMAX = 21; %% Kuhn-Tucker multiplier upper angle difference limit (u/degree)
function [ref, pv, pq] = bustypes(bus, gen)
%BUSTYPES Builds index lists for each type of bus (REF, PV, PQ).
% [REF, PV, PQ] = BUSTYPES(BUS, GEN)
% Generators with "out-of-service" status are treated as PQ buses with
% zero generation (regardless of Pg/Qg values in gen). Expects BUS and
% GEN have been converted to use internal consecutive bus numbering.
% MATPOWER
% Copyright (c) 1996-2016, Power Systems 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.
%% 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;
[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;
%% get generator status
% bus_gen_status = zeros(size(bus, 1), 1);
% bus_gen_status(gen(:, GEN_BUS)) = gen(:, GEN_STATUS) > 0;
nb = size(bus, 1);
ng = size(gen, 1);
Cg = sparse(gen(:, GEN_BUS), (1:ng)', gen(:, GEN_STATUS) > 0, nb, ng); %% gen connection matrix
%% element i, j is 1 if, generator j at bus i is ON
bus_gen_status = Cg * ones(ng, 1); %% number of generators at each bus that are ON
%% form index lists for slack, PV, and PQ buses
ref = find(bus(:, BUS_TYPE) == REF & bus_gen_status); %% reference bus index
pv = find(bus(:, BUS_TYPE) == PV & bus_gen_status); %% PV bus indices
pq = find(bus(:, BUS_TYPE) == PQ | ~bus_gen_status); %% PQ bus indices
%% pick a new reference bus if for some reason there is none (may have been shut down)
if isempty(ref)
if isempty(pv)
%% no PV bus left to convert to reference bus
else
ref = pv(1); %% use the first PV bus
pv(1) = []; %% delete it from PV list
end
end
function [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch)
%MAKEYBUS Builds the bus admittance matrix and branch admittance matrices.
% [YBUS, YF, YT] = MAKEYBUS(MPC)
% [YBUS, YF, YT] = MAKEYBUS(BASEMVA, BUS, BRANCH)
%
% Returns the full bus admittance matrix (i.e. for all buses) and the
% matrices YF and YT which, when multiplied by a complex voltage vector,
% yield the vector currents injected into each line from the "from" and
% "to" buses respectively of each line. Does appropriate conversions to p.u.
% Inputs can be a MATPOWER case struct or individual BASEMVA, BUS and
% BRANCH values. Bus numbers must be consecutive beginning at 1
% (i.e. internal ordering).
%
% See also MAKEJAC, MAKESBUS, EXT2INT.
% MATPOWER
% Copyright (c) 1996-2016, Power Systems 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.
%% extract from MPC if necessary
if nargin < 3
mpc = baseMVA;
baseMVA = mpc.baseMVA;
bus = mpc.bus;
branch = mpc.branch;
end
%% constants
nb = size(bus, 1); %% number of buses
nl = size(branch, 1); %% number of lines
%% define named indices into bus, 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;
%% check that bus numbers are equal to indices to bus (one set of bus numbers)
if any(bus(:, BUS_I) ~= (1:nb)')
error('makeYbus: buses must be numbered consecutively in bus matrix; use ext2int() to convert to internal ordering')
end
%% for each branch, compute the elements of the branch admittance matrix where
%%
%% | If | | Yff Yft | | Vf |
%% | | = | | * | |
%% | It | | Ytf Ytt | | Vt |
%%
stat = branch(:, BR_STATUS); %% ones at in-service branches
Ys = stat ./ (branch(:, BR_R) + 1j * branch(:, BR_X)); %% series admittance
Bc = stat .* branch(:, BR_B); %% line charging susceptance
tap = ones(nl, 1); %% default tap ratio = 1
i = find(branch(:, TAP)); %% indices of non-zero tap ratios
tap(i) = branch(i, TAP); %% assign non-zero tap ratios
tap = tap .* exp(1j*pi/180 * branch(:, SHIFT)); %% add phase shifters
Ytt = Ys + 1j*Bc/2;
Yff = Ytt ./ (tap .* conj(tap));
Yft = - Ys ./ conj(tap);
Ytf = - Ys ./ tap;
%% compute shunt admittance
%% if Psh is the real power consumed by the shunt at V = 1.0 p.u.
%% and Qsh is the reactive power injected by the shunt at V = 1.0 p.u.
%% then Psh - j Qsh = V * conj(Ysh * V) = conj(Ysh) = Gs - j Bs,
%% i.e. Ysh = Psh + j Qsh, so ...
Ysh = (bus(:, GS) + 1j * bus(:, BS)) / baseMVA; %% vector of shunt admittances
%% build connection matrices
f = branch(:, F_BUS); %% list of "from" buses
t = branch(:, T_BUS); %% list of "to" buses
Cf = sparse(1:nl, f, ones(nl, 1), nl, nb); %% connection matrix for line & from buses
Ct = sparse(1:nl, t, ones(nl, 1), nl, nb); %% connection matrix for line & to buses
%% build Yf and Yt such that Yf * V is the vector of complex branch currents injected
%% at each branch's "from" bus, and Yt is the same for the "to" bus end
i = [1:nl; 1:nl]'; %% double set of row indices
Yf = sparse(i, [f; t], [Yff; Yft], nl, nb);
Yt = sparse(i, [f; t], [Ytf; Ytt], nl, nb);
% Yf = spdiags(Yff, 0, nl, nl) * Cf + spdiags(Yft, 0, nl, nl) * Ct;
% Yt = spdiags(Ytf, 0, nl, nl) * Cf + spdiags(Ytt, 0, nl, nl) * Ct;
%% build Ybus
Ybus = Cf' * Yf + Ct' * Yt + ... %% branch admittances
sparse(1:nb, 1:nb, Ysh, nb, nb); %% shunt admittance
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
% p.u., and overrides the PG and QG columns in GEN, using GEN only for
% connectivity information.
%
% [SBUS, DSBUS_DVM] = MAKESBUS(BASEMVA, BUS, GEN, MPOPT, VM)
% With two output arguments, it computes the partial derivative of the
% bus injections with respect to voltage magnitude, leaving the first
% return value SBUS empty. If VM is empty, it assumes no voltage dependence
% and returns a sparse zero matrix.
%
% See also MAKEYBUS.
% MATPOWER
% Copyright (c) 1996-2016, Power Systems 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 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 [V, converged, i,J] = 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, Power Systems 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.
mpopt.pf.nr.max_it=20;
%% 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));