|
| 1 | +function t_opf_dc_knitro(quiet) |
| 2 | +% t_opf_dc_knitro - Tests for legacy DC optimal power flow using Artelys Knitro solver. |
| 3 | + |
| 4 | +% MATPOWER |
| 5 | +% Copyright (c) 2004-2025, Power Systems Engineering Research Center (PSERC) |
| 6 | +% by Ray Zimmerman, PSERC Cornell |
| 7 | +% |
| 8 | +% This file is part of MATPOWER. |
| 9 | +% Covered by the 3-clause BSD License (see LICENSE file for details). |
| 10 | +% See https://matpower.org for more info. |
| 11 | + |
| 12 | +if nargin < 1 |
| 13 | + quiet = 0; |
| 14 | +end |
| 15 | + |
| 16 | +num_tests = 43; |
| 17 | + |
| 18 | +t_begin(num_tests, quiet); |
| 19 | + |
| 20 | +[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... |
| 21 | + VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; |
| 22 | +[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... |
| 23 | + MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... |
| 24 | + QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; |
| 25 | +[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... |
| 26 | + TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... |
| 27 | + ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; |
| 28 | + |
| 29 | +casefile = 't_case9_opf'; |
| 30 | +if quiet |
| 31 | + verbose = 0; |
| 32 | +else |
| 33 | + verbose = 0; |
| 34 | +end |
| 35 | +if have_feature('octave') |
| 36 | + if have_feature('octave', 'vnum') >= 4 |
| 37 | + file_in_path_warn_id = 'Octave:data-file-in-path'; |
| 38 | + else |
| 39 | + file_in_path_warn_id = 'Octave:load-file-in-path'; |
| 40 | + end |
| 41 | + s1 = warning('query', file_in_path_warn_id); |
| 42 | + warning('off', file_in_path_warn_id); |
| 43 | +end |
| 44 | + |
| 45 | +t0 = 'DC OPF (Knitro): '; |
| 46 | +mpopt = mpoption('out.all', 0, 'verbose', verbose); |
| 47 | +mpopt = mpoption(mpopt, 'opf.dc.solver', 'KNITRO'); |
| 48 | + |
| 49 | +if have_feature('knitro') |
| 50 | + mpopt = mpoption(mpopt, 'knitro.opts.opttol', 1e-10); |
| 51 | + % mpopt = mpoption(mpopt, 'knitro.opts.opttol_abs', 1e-8); |
| 52 | + |
| 53 | + %% set up indices |
| 54 | + ib_data = [1:BUS_AREA BASE_KV:VMIN]; |
| 55 | + ib_voltage = [VM VA]; |
| 56 | + ib_lam = [LAM_P LAM_Q]; |
| 57 | + ib_mu = [MU_VMAX MU_VMIN]; |
| 58 | + ig_data = [GEN_BUS QMAX QMIN MBASE:APF]; |
| 59 | + ig_disp = [PG QG VG]; |
| 60 | + ig_mu = (MU_PMAX:MU_QMIN); |
| 61 | + ibr_data = (1:ANGMAX); |
| 62 | + ibr_flow = (PF:QT); |
| 63 | + ibr_mu = [MU_SF MU_ST]; |
| 64 | + ibr_angmu = [MU_ANGMIN MU_ANGMAX]; |
| 65 | + |
| 66 | + %% get solved DC power flow case from MAT-file |
| 67 | + load soln9_dcopf; %% defines bus_soln, gen_soln, branch_soln, f_soln |
| 68 | + |
| 69 | + %% run OPF |
| 70 | + t = t0; |
| 71 | + [baseMVA, bus, gen, gencost, branch, f, success, et] = rundcopf(casefile, mpopt); |
| 72 | + t_ok(success, [t 'success']); |
| 73 | + t_is(f, f_soln, 3, [t 'f']); |
| 74 | + t_is( bus(:,ib_data ), bus_soln(:,ib_data ), 10, [t 'bus data']); |
| 75 | + t_is( bus(:,ib_voltage), bus_soln(:,ib_voltage), 3, [t 'bus voltage']); |
| 76 | + t_is( bus(:,ib_lam ), bus_soln(:,ib_lam ), 3, [t 'bus lambda']); |
| 77 | + t_is( bus(:,ib_mu ), bus_soln(:,ib_mu ), 2, [t 'bus mu']); |
| 78 | + t_is( gen(:,ig_data ), gen_soln(:,ig_data ), 10, [t 'gen data']); |
| 79 | + t_is( gen(:,ig_disp ), gen_soln(:,ig_disp ), 3, [t 'gen dispatch']); |
| 80 | + t_is( gen(:,ig_mu ), gen_soln(:,ig_mu ), 3, [t 'gen mu']); |
| 81 | + t_is(branch(:,ibr_data ), branch_soln(:,ibr_data ), 10, [t 'branch data']); |
| 82 | + t_is(branch(:,ibr_flow ), branch_soln(:,ibr_flow ), 3, [t 'branch flow']); |
| 83 | + t_is(branch(:,ibr_mu ), branch_soln(:,ibr_mu ), 2, [t 'branch mu']); |
| 84 | + |
| 85 | + %%----- test OPF with angle difference limits ----- |
| 86 | + t = [t0 'w/angle diff lims : ']; |
| 87 | + mpc = loadcase(casefile); |
| 88 | + mpc.branch(4, ANGMAX) = 3; |
| 89 | + mpc.branch(7, ANGMIN) = -4.5; |
| 90 | + r = rundcopf(mpc, mpopt); |
| 91 | + [bus, gen, branch, f, success] = deal(r.bus, r.gen, r.branch, r.f, r.success); |
| 92 | + t_ok(success, [t 'success']); |
| 93 | + t_is( f, 6456.7213, 3, [t 'f']); |
| 94 | + t_is( bus(:,ib_data ), bus_soln(:,ib_data ), 10, [t 'bus data']); |
| 95 | + t_is( gen(:,ig_data ), gen_soln(:,ig_data ), 10, [t 'gen data']); |
| 96 | + t_is( gen(:,PG ), [99.98497;89.35133;125.66371], 4, [t 'gen dispatch']); |
| 97 | + t_is(branch(:,ibr_data ), mpc.branch(:,ibr_data ), 10, [t 'branch data']); |
| 98 | + e = zeros(size(branch, 1), 1); |
| 99 | + e(4) = 297.83776; |
| 100 | + e(7) = -26.94788; |
| 101 | + t_is(branch(:,MU_ANGMAX )-branch(:,MU_ANGMIN ), e, 4, [t 'branch ang diff mu']); |
| 102 | + |
| 103 | + t = [t0 'w/ignored angle diff lims : ']; |
| 104 | + mpopt1 = mpoption(mpopt, 'opf.ignore_angle_lim', 1); |
| 105 | + r = rundcopf(mpc, mpopt1); |
| 106 | + [bus, gen, branch, f, success] = deal(r.bus, r.gen, r.branch, r.f, r.success); |
| 107 | + t_ok(success, [t 'success']); |
| 108 | + t_is(f, f_soln, 3, [t 'f']); |
| 109 | + t_is( bus(:,ib_data ), bus_soln(:,ib_data ), 10, [t 'bus data']); |
| 110 | + t_is( bus(:,ib_voltage), bus_soln(:,ib_voltage), 3, [t 'bus voltage']); |
| 111 | + t_is( bus(:,ib_lam ), bus_soln(:,ib_lam ), 3, [t 'bus lambda']); |
| 112 | + t_is( bus(:,ib_mu ), bus_soln(:,ib_mu ), 2, [t 'bus mu']); |
| 113 | + t_is( gen(:,ig_data ), gen_soln(:,ig_data ), 10, [t 'gen data']); |
| 114 | + t_is( gen(:,ig_disp ), gen_soln(:,ig_disp ), 3, [t 'gen dispatch']); |
| 115 | + t_is( gen(:,ig_mu ), gen_soln(:,ig_mu ), 3, [t 'gen mu']); |
| 116 | + t_is(branch(:,ibr_data ), mpc.branch(:,ibr_data ), 10, [t 'branch data']); |
| 117 | + t_is(branch(:,ibr_flow ), branch_soln(:,ibr_flow ), 3, [t 'branch flow']); |
| 118 | + t_is(branch(:,ibr_mu ), branch_soln(:,ibr_mu ), 2, [t 'branch mu']); |
| 119 | + |
| 120 | + %%----- run OPF with extra linear user constraints & costs ----- |
| 121 | + %% two new z variables |
| 122 | + %% 0 <= z1, P3 - P1 <= z1 |
| 123 | + %% 0 <= z2, P3 - P2 <= z2 |
| 124 | + %% with A and N sized for DC opf |
| 125 | + mpc = loadcase(casefile); |
| 126 | + mpc.A = sparse([1;1;1;2;2;2],[10;12;13;12;11;14],[-1;1;-1;1;-1;-1],2,14); |
| 127 | + mpc.u = [0; 0]; |
| 128 | + mpc.l = [-Inf; -Inf]; |
| 129 | + mpc.zl = [0; 0]; |
| 130 | + |
| 131 | + mpc.N = sparse([1;2], [13;14], [1;1], 2, 14); %% new z variables only |
| 132 | + mpc.fparm = ones(2,1) * [1 0 0 1]; %% w = r = z |
| 133 | + mpc.H = sparse(2,2); %% no quadratic term |
| 134 | + mpc.Cw = [1000;1]; |
| 135 | + |
| 136 | + t = [t0 'w/extra constraints & costs 1 : ']; |
| 137 | + [r, success] = rundcopf(mpc, mpopt); |
| 138 | + t_ok(success, [t 'success']); |
| 139 | + t_is(r.gen(1, PG), 116.15974, 4, [t 'Pg1 = 116.15974']); |
| 140 | + t_is(r.gen(3, PG), 116.15974, 4, [t 'Pg3 = 116.15974']); |
| 141 | + t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']); |
| 142 | + t_is(r.cost.usr, 0.3348, 3, [t 'user costs']); |
| 143 | + |
| 144 | + %% with A and N sized for AC opf |
| 145 | + mpc = loadcase(casefile); |
| 146 | + mpc.A = sparse([1;1;1;2;2;2],[19;21;25;21;20;26],[-1;1;-1;1;-1;-1],2,26); |
| 147 | + mpc.u = [0; 0]; |
| 148 | + mpc.l = [-Inf; -Inf]; |
| 149 | + mpc.zl = [0; 0]; |
| 150 | + |
| 151 | + mpc.N = sparse([1;2], [25;26], [1;1], 2, 26); %% new z variables only |
| 152 | + mpc.fparm = ones(2,1) * [1 0 0 1]; %% w = r = z |
| 153 | + mpc.H = sparse(2,2); %% no quadratic term |
| 154 | + mpc.Cw = [1000;1]; |
| 155 | + |
| 156 | + t = [t0 'w/extra constraints & costs 2 : ']; |
| 157 | + [r, success] = rundcopf(mpc, mpopt); |
| 158 | + t_ok(success, [t 'success']); |
| 159 | + t_is(r.gen(1, PG), 116.15974, 4, [t 'Pg1 = 116.15974']); |
| 160 | + t_is(r.gen(3, PG), 116.15974, 4, [t 'Pg3 = 116.15974']); |
| 161 | + t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']); |
| 162 | + t_is(r.cost.usr, 0.3348, 3, [t 'user costs']); |
| 163 | + |
| 164 | + t = [t0 'infeasible : ']; |
| 165 | + %% with A and N sized for DC opf |
| 166 | + mpc = loadcase(casefile); |
| 167 | + mpc.A = sparse([1;1], [10;11], [1;1], 1, 14); %% Pg1 + Pg2 |
| 168 | + mpc.u = Inf; |
| 169 | + mpc.l = 600; |
| 170 | + [r, success] = rundcopf(mpc, mpopt); |
| 171 | + t_ok(~success, [t 'no success']); |
| 172 | + |
| 173 | + %% OPF with all buses isolated |
| 174 | + t = [t0 'all buses isolated : ']; |
| 175 | + mpc = loadcase(casefile); |
| 176 | + mpc.bus(:, BUS_TYPE) = NONE; |
| 177 | + try |
| 178 | + r = rundcopf(mpc, mpopt); |
| 179 | + t_is(r.success, 0, 12, [t 'success = 0']); |
| 180 | + catch |
| 181 | + t_ok(0, [t 'unexpected fatal error']); |
| 182 | + end |
| 183 | +else |
| 184 | + t_skip(num_tests, 'Artelys Knitro not available'); |
| 185 | +end |
| 186 | + |
| 187 | +if have_feature('octave') |
| 188 | + warning(s1.state, file_in_path_warn_id); |
| 189 | +end |
| 190 | + |
| 191 | +t_end; |
0 commit comments