0001 function t_opf_dc_glpk(quiet)
0002
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 if nargin < 1
0034 quiet = 0;
0035 end
0036
0037 algs = [1; 1; 2];
0038 if have_fcn('octave')
0039 v = ver('Octave');
0040 if vstr2num(v.Version) < 3.007
0041 dual = [0; 1; 1];
0042 else
0043 dual = [1; 2; 1];
0044 end
0045 else
0046 dual = [0; 2; 1];
0047 end
0048 num_tests = 23 * length(algs);
0049
0050 t_begin(num_tests, quiet);
0051
0052 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0053 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0054 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0055 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0056 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0057 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0058 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0059 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0060
0061 casefile = 't_case9_opf';
0062 if quiet
0063 verbose = 0;
0064 else
0065 verbose = 0;
0066 end
0067 if have_fcn('octave')
0068 s1 = warning('query', 'Octave:load-file-in-path');
0069 warning('off', 'Octave:load-file-in-path');
0070 end
0071
0072 s2 = warning('query', 'MATLAB:singularMatrix');
0073
0074 mpopt = mpoption('out.all', 0, 'verbose', verbose);
0075 mpopt = mpoption(mpopt, 'opf.dc.solver', 'GLPK');
0076
0077
0078 if have_fcn('glpk')
0079 for k = 1:length(algs)
0080 mpopt = mpoption(mpopt, 'glpk.opts.lpsolver', algs(k), 'glpk.opts.dual', dual(k));
0081 methods = {
0082 'primal simplex',
0083 'dual simplex',
0084 'interior',
0085 };
0086 t0 = sprintf('DC OPF (GLPK %s): ', methods{k});
0087
0088
0089 ib_data = [1:BUS_AREA BASE_KV:VMIN];
0090 ib_voltage = [VM VA];
0091 ib_lam = [LAM_P LAM_Q];
0092 ib_mu = [MU_VMAX MU_VMIN];
0093 ig_data = [GEN_BUS QMAX QMIN MBASE:APF];
0094 ig_disp = [PG QG VG];
0095 ig_mu = (MU_PMAX:MU_QMIN);
0096 ibr_data = (1:ANGMAX);
0097 ibr_flow = (PF:QT);
0098 ibr_mu = [MU_SF MU_ST];
0099 ibr_angmu = [MU_ANGMIN MU_ANGMAX];
0100
0101
0102 load soln9_dcopf;
0103
0104
0105 t = t0;
0106 [baseMVA, bus, gen, gencost, branch, f, success, et] = rundcopf(casefile, mpopt);
0107 t_ok(success, [t 'success']);
0108 t_is(f, f_soln, 3, [t 'f']);
0109 t_is( bus(:,ib_data ), bus_soln(:,ib_data ), 10, [t 'bus data']);
0110 t_is( bus(:,ib_voltage), bus_soln(:,ib_voltage), 3, [t 'bus voltage']);
0111 t_is( bus(:,ib_lam ), bus_soln(:,ib_lam ), 3, [t 'bus lambda']);
0112 t_is( bus(:,ib_mu ), bus_soln(:,ib_mu ), 2, [t 'bus mu']);
0113 t_is( gen(:,ig_data ), gen_soln(:,ig_data ), 10, [t 'gen data']);
0114 t_is( gen(:,ig_disp ), gen_soln(:,ig_disp ), 3, [t 'gen dispatch']);
0115 t_is( gen(:,ig_mu ), gen_soln(:,ig_mu ), 3, [t 'gen mu']);
0116 t_is(branch(:,ibr_data ), branch_soln(:,ibr_data ), 10, [t 'branch data']);
0117 t_is(branch(:,ibr_flow ), branch_soln(:,ibr_flow ), 3, [t 'branch flow']);
0118 t_is(branch(:,ibr_mu ), branch_soln(:,ibr_mu ), 2, [t 'branch mu']);
0119
0120
0121
0122
0123
0124
0125 mpc = loadcase(casefile);
0126 mpc.A = sparse([1;1;1;2;2;2],[10;11;13;11;12;14],[-1;1;-1;1;-1;-1],2,14);
0127 mpc.u = [0; 0];
0128 mpc.l = [-Inf; -Inf];
0129 mpc.zl = [0; 0];
0130
0131 mpc.N = sparse([1;2], [13;14], [1;1], 2, 14);
0132 mpc.fparm = ones(2,1) * [1 0 0 1];
0133 mpc.H = sparse(2,2);
0134 mpc.Cw = [1000;1];
0135
0136 t = [t0 'w/extra constraints & costs 1 : '];
0137 [r, success] = rundcopf(mpc, mpopt);
0138 t_ok(success, [t 'success']);
0139 t_is(r.gen(1, PG), 116.15974, 5, [t 'Pg1 = 116.15974']);
0140 t_is(r.gen(2, PG), 116.15974, 5, [t 'Pg2 = 116.15974']);
0141 t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']);
0142 t_is(r.cost.usr, 0.3348, 4, [t 'user costs']);
0143
0144
0145 mpc = loadcase(casefile);
0146 mpc.A = sparse([1;1;1;2;2;2],[19;20;25;20;21;26],[-1;1;-1;1;-1;-1],2,26);
0147 mpc.u = [0; 0];
0148 mpc.l = [-Inf; -Inf];
0149 mpc.zl = [0; 0];
0150
0151 mpc.N = sparse([1;2], [25;26], [1;1], 2, 26);
0152 mpc.fparm = ones(2,1) * [1 0 0 1];
0153 mpc.H = sparse(2,2);
0154 mpc.Cw = [1000;1];
0155
0156 t = [t0 'w/extra constraints & costs 2 : '];
0157 [r, success] = rundcopf(mpc, mpopt);
0158 t_ok(success, [t 'success']);
0159 t_is(r.gen(1, PG), 116.15974, 5, [t 'Pg1 = 116.15974']);
0160 t_is(r.gen(2, PG), 116.15974, 5, [t 'Pg2 = 116.15974']);
0161 t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']);
0162 t_is(r.cost.usr, 0.3348, 4, [t 'user costs']);
0163
0164 t = [t0 'infeasible : '];
0165 warning('off', 'MATLAB:singularMatrix');
0166
0167 mpc = loadcase(casefile);
0168 mpc.A = sparse([1;1], [10;11], [1;1], 1, 14);
0169 mpc.u = Inf;
0170 mpc.l = 600;
0171 [r, success] = rundcopf(mpc, mpopt);
0172 t_ok(~success, [t 'no success']);
0173
0174 end
0175 else
0176 t_skip(num_tests, 'GLPK not available');
0177 end
0178
0179 if have_fcn('octave')
0180 warning(s1.state, 'Octave:load-file-in-path');
0181 end
0182 warning(s2.state, 'MATLAB:singularMatrix');
0183
0184 t_end;
0185
0186
0187 function num = vstr2num(vstr)
0188
0189
0190 pat = '\.?(\d+)';
0191 [s,e,tE,m,t] = regexp(vstr, pat);
0192 b = 1;
0193 num = 0;
0194 for k = 1:length(t)
0195 num = num + b * str2num(t{k}{1});
0196 b = b / 1000;
0197 end