t_vdep_load

T_VDEP_LOAD    Test voltage dependent ZIP load model for PF, CPF, OPF.

0001 function t_vdep_load(quiet)
0002 %T_VDEP_LOAD    Test voltage dependent ZIP load model for PF, CPF, OPF.
0003 
0004 %   MATPOWER
0005 %   Copyright (c) 2009-2016, Power Systems Engineering Research Center (PSERC)
0006 %   by Ray Zimmerman, PSERC Cornell
0007 %
0008 %   This file is part of MATPOWER.
0009 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0010 %   See https://matpower.org for more info.
0011 
0012 if nargin < 1
0013     quiet = 0;
0014 end
0015 
0016 pfalgs = {'NR', 'FDXB', 'FDBX'};
0017 cpfparm = {'natural', 'arc-len', 'pseudo arc-len'};
0018 
0019 t_begin(33+27*length(pfalgs)+39*2*length(cpfparm), quiet);
0020 
0021 if quiet
0022     verbose = 0;
0023 else
0024     verbose = 0;
0025 end
0026 verbose = 0;
0027 
0028 casefile = 't_case9_opf';
0029 mpopt = mpoption;
0030 mpopt = mpoption(mpopt, 'out.gen', 1);
0031 mpopt = mpoption(mpopt, 'pf.tol', 1e-9);
0032 mpopt = mpoption(mpopt, 'opf.ac.solver', 'MIPS');
0033 mpopt = mpoption(mpopt, 'opf.violation', 1e-6, 'mips.gradtol', 1e-9, ...
0034         'mips.comptol', 1e-8, 'mips.costtol', 1e-8);
0035 mpopt = mpoption(mpopt, 'verbose', verbose);
0036 if ~verbose
0037     mpopt = mpoption(mpopt, 'out.all', 0);
0038 end
0039 
0040 %% define named indices into data matrices
0041 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0042     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0043 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0044     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0045     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0046 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0047     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0048     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0049 
0050 %% load base case file
0051 mpc = loadcase(casefile);
0052 nb = size(mpc.bus, 1);
0053 
0054 t = 'PF - base case : ';
0055 mpopt0 = mpopt;
0056 r0 = runpf(mpc, mpopt);
0057 t_ok(r0.success, [t 'success']);
0058 
0059 for k = 1:length(pfalgs)
0060     mpopt = mpoption(mpopt0, 'pf.alg', pfalgs{k});
0061 
0062     t = sprintf('PF (%s) - constant power only : ', pfalgs{k});
0063     zip_w = [1 0 0];
0064     mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads', struct('pw', zip_w, 'qw', zip_w));
0065     r = runpf(mpc, mpopt);
0066     t_ok(r.success, [t 'success']);
0067     t_is(r.bus, r0.bus, 6, [t 'bus']);
0068     t_is(r.gen, r0.gen, 6, [t 'gen']);
0069     t_is(r.branch, r0.branch, 6, [t 'branch']);
0070 
0071     t = sprintf('PF (%s) - constant current only : ', pfalgs{k});
0072     zip_w = [0 1 0];
0073     mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads', struct('pw', zip_w, 'qw', zip_w));
0074     r = runpf(mpc, mpopt);
0075     [tPd, tQd] = total_load(r, 'bus', [], mpopt);
0076     % [tPd, tQd] = total_load(r, 'bus');
0077     mpc1 = mpc;
0078     mpc1.bus(:, PD) = mpc1.bus(:, PD) .* r.bus(:, VM);
0079     mpc1.bus(:, QD) = mpc1.bus(:, QD) .* r.bus(:, VM);
0080     r1 = runpf(mpc1, mpopt0);
0081     t_ok(r1.success, [t 'success']);
0082     t_ok(r.success, [t 'success']);
0083     t_is(r.bus(:, [VM VA]), r1.bus(:, [VM VA]), 6, [t 'bus voltages']);
0084     t_is(tPd, r1.bus(:, PD), 6, [t 'bus P loads']);
0085     t_is(tPd, [0; 0; 0; 0; 87.937521; 0; 98.674071; 0; 120.041555], 6, [t 'bus P loads']);
0086     t_is(tQd, r1.bus(:, QD), 6, [t 'bus Q loads']);
0087     t_is(tQd, [0; 0; 0; 0; 29.312507; 0; 34.535925; 0; 48.016622], 6, [t 'bus Q loads']);
0088     t_is(r.gen(:, [PG QG]), r1.gen(:, [PG QG]), 6, [t 'gen dispatches']);
0089     t_is(r.branch, r1.branch, 6, [t 'branch']);
0090 
0091     t = sprintf('PF (%s) - constant impedance only : ', pfalgs{k});
0092     mpc1 = mpc;
0093     mpc1.bus(:, GS) =  mpc1.bus(:, PD);
0094     mpc1.bus(:, BS) = -mpc1.bus(:, QD);
0095     mpc1.bus(:, PD) = 0;
0096     mpc1.bus(:, QD) = 0;
0097     r1 = runpf(mpc1, mpopt0);
0098     t_ok(r1.success, [t 'success']);
0099 
0100     zip_w = [0 0 1];
0101     mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads', struct('pw', zip_w, 'qw', []));
0102     r = runpf(mpc, mpopt);
0103     t_ok(r.success, [t 'success']);
0104     t_is(r.bus(:, [VM VA]), r1.bus(:, [VM VA]), 6, [t 'bus voltages']);
0105     t_is(r.gen(:, [PG QG]), r1.gen(:, [PG QG]), 6, [t 'gen dispatches']);
0106     t_is(r.branch, r1.branch, 6, [t 'branch']);
0107 
0108     t = sprintf('PF (%s) - combo ZIP loads : ', pfalgs{k});
0109     zip_w = [0.1 0.4 0.5];
0110     mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads.pw', zip_w);
0111     r = runpf(mpc, mpopt);
0112     [tPd, tQd] = total_load(r, 'bus', [], mpopt);
0113     mpc1 = mpc;
0114     Vm = r.bus(:, VM);
0115     scale = [Vm.^0 Vm Vm.^2] * zip_w';
0116     mpc1.bus(:, PD) = mpc1.bus(:, PD) .* scale;
0117     mpc1.bus(:, QD) = mpc1.bus(:, QD) .* scale;
0118     r1 = runpf(mpc1, mpopt0);
0119     t_ok(r1.success, [t 'success']);
0120 
0121     t_ok(r.success, [t 'success']);
0122     t_is(r.bus(:, [VM VA]), r1.bus(:, [VM VA]), 6, [t 'bus voltages']);
0123     t_is(tPd, r1.bus(:, PD), 6, [t 'bus P loads']);
0124     t_is(tPd, [0; 0; 0; 0; 87.205063; 0; 98.205247; 0; 118.314510], 6, [t 'bus P loads']);
0125     t_is(tQd, r1.bus(:, QD), 6, [t 'bus Q loads']);
0126     t_is(tQd, [0; 0; 0; 0; 29.068354; 0; 34.371836; 0; 47.325804], 6, [t 'bus Q loads']);
0127     t_is(r.gen(:, [PG QG]), r1.gen(:, [PG QG]), 6, [t 'gen dispatches']);
0128     t_is(r.branch, r1.branch, 6, [t 'branch']);
0129 end
0130 
0131 t = 'OPF - base case : ';
0132 % mpopt0 = mpoption(mpopt0, 'verbose', 2);
0133 mpopt = mpopt0;
0134 r0 = runopf(mpc, mpopt);
0135 t_ok(r0.success, [t 'success']);
0136 
0137 t = 'OPF - constant power only : ';
0138 zip_w = [1 0 0];
0139 mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads.pw', zip_w);
0140 r = runopf(mpc, mpopt);
0141 t_ok(r.success, [t 'success']);
0142 t_is(r.f, r0.f, 6, [t 'f']);
0143 t_is(r.bus, r0.bus, 6, [t 'bus']);
0144 t_is(r.gen, r0.gen, 6, [t 'gen']);
0145 t_is(r.branch, r0.branch, 6, [t 'branch']);
0146 
0147 t = 'OPF - constant current only : ';
0148 zip_w = [0 1 0];
0149 mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads.pw', zip_w);
0150 r = runopf(mpc, mpopt);
0151 [tPd, tQd] = total_load(r, 'bus', [], mpopt);
0152 % [tPd, tQd] = total_load(r, 'bus');
0153 mpc1 = mpc;
0154 mpc1.bus(:, PD) = mpc1.bus(:, PD) .* r.bus(:, VM);
0155 mpc1.bus(:, QD) = mpc1.bus(:, QD) .* r.bus(:, VM);
0156 k = find(mpc1.bus(:, PD));  %% buses with non-zero loads
0157 mpc1.bus(k, VMAX) = r.bus(k, VM);
0158 mpc1.bus(k, VMIN) = r.bus(k, VM);
0159 r1 = runopf(mpc1, mpopt0);
0160 t_ok(r1.success, [t 'success']);
0161 t_ok(r.success, [t 'success']);
0162 t_is(r.f, r1.f, 4, [t 'f']);
0163 t_is(r.bus(:, [VM VA]), r1.bus(:, [VM VA]), 6, [t 'bus voltages']);
0164 t_is(tPd, r1.bus(:, PD), 6, [t 'bus P loads']);
0165 t_is(tPd, [0; 0; 0; 0; 81.708978; 0; 90.914185; 0; 112.500000], 6, [t 'bus P loads']);
0166 t_is(tQd, r1.bus(:, QD), 6, [t 'bus Q loads']);
0167 t_is(tQd, [0; 0; 0; 0; 27.236326; 0; 31.819965; 0; 45.000000], 6, [t 'bus Q loads']);
0168 t_is(r.gen(:, [PG QG]), r1.gen(:, [PG QG]), 6, [t 'gen dispatches']);
0169 t_is(r.branch(:, 1:ANGMAX), r1.branch(:, 1:ANGMAX), 6, [t 'branch']);
0170 
0171 t = 'OPF - constant impedance only : ';
0172 mpc1 = mpc;
0173 mpc1.bus(:, GS) =  mpc1.bus(:, PD);
0174 mpc1.bus(:, BS) = -mpc1.bus(:, QD);
0175 mpc1.bus(:, PD) = 0;
0176 mpc1.bus(:, QD) = 0;
0177 r1 = runopf(mpc1, mpopt0);
0178 t_ok(r1.success, [t 'success']);
0179 
0180 zip_w = [0 0 1];
0181 mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads.pw', zip_w);
0182 r = runopf(mpc, mpopt);
0183 t_ok(r.success, [t 'success']);
0184 t_is(r.f, r1.f, 6, [t 'f']);
0185 t_is(r.bus(:, [VM VA]), r1.bus(:, [VM VA]), 6, [t 'bus voltages']);
0186 t_is(r.gen(:, [PG QG]), r1.gen(:, [PG QG]), 6, [t 'gen dispatches']);
0187 t_is(r.branch, r1.branch, 6, [t 'branch']);
0188 
0189 t = 'OPF - combo ZIP loads : ';
0190 zip_w = [0.1 0.4 0.5];
0191 mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads.pw', zip_w);
0192 r = runopf(mpc, mpopt);
0193 [tPd, tQd] = total_load(r, 'bus', [], mpopt);
0194 mpc1 = mpc;
0195 Vm = r.bus(:, VM);
0196 scale = [Vm.^0 Vm Vm.^2] * zip_w';
0197 mpc1.bus(:, PD) = mpc1.bus(:, PD) .* scale;
0198 mpc1.bus(:, QD) = mpc1.bus(:, QD) .* scale;
0199 k = find(mpc1.bus(:, PD));  %% buses with non-zero loads
0200 mpc1.bus(k, VMAX) = r.bus(k, VM);
0201 mpc1.bus(k, VMIN) = r.bus(k, VM);
0202 r1 = runopf(mpc1, mpopt0);
0203 t_ok(r1.success, [t 'success']);
0204 
0205 t_ok(r.success, [t 'success']);
0206 t_is(r.f, r1.f, 4, [t 'f']);
0207 t_is(r.bus(:, [VM VA]), r1.bus(:, [VM VA]), 6, [t 'bus voltages']);
0208 t_is(tPd, r1.bus(:, PD), 6, [t 'bus P loads']);
0209 t_is(tPd, [0; 0; 0; 0; 78.909380; 0; 87.300375; 0; 108.125000], 6, [t 'bus P loads']);
0210 t_is(tQd, r1.bus(:, QD), 6, [t 'bus Q loads']);
0211 t_is(tQd, [0; 0; 0; 0; 26.303127; 0; 30.555131; 0; 43.250000], 6, [t 'bus Q loads']);
0212 t_is(r.gen(:, [PG QG]), r1.gen(:, [PG QG]), 6, [t 'gen dispatches']);
0213 t_is(r.branch(:, 1:ANGMAX), r1.branch(:, 1:ANGMAX), 6, [t 'branch']);
0214 
0215 
0216 mpopt =  mpopt0;
0217     plot_nose_curve = 0;
0218     mpopt = mpoption(mpopt, 'cpf.stop_at', 1);
0219     mpopt = mpoption(mpopt, 'cpf.step', 0.05);
0220     %mpopt = mpoption(mpopt, 'cpf.adapt_step', 1);
0221     %mpopt = mpoption(mpopt, 'cpf.adapt_step_tol', 2e-5);
0222     mpopt = mpoption(mpopt, 'cpf.plot.level', plot_nose_curve);
0223     %mpopt = mpoption(mpopt, 'cpf.plot.bus', 9);
0224 mpopt0 = mpopt;
0225 
0226     %% set up base and target cases
0227     mpcb = loadcase(casefile);
0228     mpct = mpcb;
0229     factor = 2.5 * 0.7;
0230     mpct.gen(:, [PG QG]) = mpct.gen(:, [PG QG]) * factor;
0231     mpct.bus(:, [PD QD]) = mpct.bus(:, [PD QD]) * factor;
0232 
0233     iterations = [
0234         20 20 20 20;
0235         21 22 22 22;
0236         21 22 22 22;
0237     ];
0238 
0239 for k = 1:length(cpfparm)
0240     mpopt0 = mpoption(mpopt0, 'cpf.parameterization', k);
0241 
0242     t = sprintf('CPF (%s) - base case : ', cpfparm{k});
0243     r0 = runcpf(mpcb, mpct, mpopt0);
0244     t_ok(r0.success, [t 'success']);
0245     t_is(r0.cpf.iterations, iterations(k, 1), 12, [t 'iterations']);
0246     t_is(r0.cpf.max_lam, 1, 12, [t 'max_lam']);
0247 
0248     mpopt = mpopt0;
0249     t = sprintf('CPF (%s) - constant power only : ', cpfparm{k});
0250     zip_w = [1 0 0];
0251     mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads.pw', zip_w);
0252     r = runcpf(mpcb, mpct, mpopt);
0253     t_ok(r.success, [t 'success']);
0254     t_is(r.cpf.iterations, r0.cpf.iterations, 12, [t 'iterations']);
0255     t_is(r.cpf.max_lam, r0.cpf.max_lam, 12, [t 'max_lam']);
0256     t_is(r.bus, r0.bus, 6, [t 'bus']);
0257     t_is(r.gen, r0.gen, 6, [t 'gen']);
0258     t_is(r.branch, r0.branch, 6, [t 'branch']);
0259 
0260     t = sprintf('CPF (%s) - constant current only : ', cpfparm{k});
0261     zip_w = [0 1 0];
0262     mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads.pw', zip_w);
0263     rb = runpf(mpcb, mpopt);
0264     r = runcpf(mpcb, mpct, mpopt);
0265     [tPd, tQd] = total_load(r, 'bus', [], mpopt);
0266     mpc1b = mpcb;
0267     mpc1b.bus(:, PD) = mpc1b.bus(:, PD) .* rb.bus(:, VM);
0268     mpc1b.bus(:, QD) = mpc1b.bus(:, QD) .* rb.bus(:, VM);
0269     mpc1t = mpct;
0270     mpc1t.bus(:, PD) = mpc1t.bus(:, PD) .* r.bus(:, VM);
0271     mpc1t.bus(:, QD) = mpc1t.bus(:, QD) .* r.bus(:, VM);
0272     r1 = runcpf(mpc1b, mpc1t, mpopt0);
0273     t_ok(r1.success, [t 'success']);
0274     t_is(r1.cpf.iterations, iterations(k,2), 12, [t 'iterations']);
0275     t_ok(r.success, [t 'success']);
0276     t_is(r.bus(:, [VM VA]), r1.bus(:, [VM VA]), 6, [t 'bus voltages']);
0277     t_is(tPd, r1.bus(:, PD), 6, [t 'bus P loads']);
0278     t_is(tPd, [0; 0; 0; 0; 143.494121; 0; 163.580859; 0; 191.911932], 6, [t 'bus P loads']);
0279     t_is(tQd, r1.bus(:, QD), 6, [t 'bus Q loads']);
0280     t_is(tQd, [0; 0; 0; 0; 47.831374; 0; 57.253301; 0; 76.764773], 6, [t 'bus Q loads']);
0281     t_is(r.gen(:, [PG QG]), r1.gen(:, [PG QG]), 6, [t 'gen dispatches']);
0282     t_is(r.branch, r1.branch, 6, [t 'branch']);
0283 
0284     t = sprintf('CPF (%s) - constant impedance only : ', cpfparm{k});
0285     zip_w = [0 0 1];
0286     mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads.pw', zip_w);
0287     rb = runpf(mpcb, mpopt);
0288     r = runcpf(mpcb, mpct, mpopt);
0289     [tPd, tQd] = total_load(r, 'bus', [], mpopt);
0290     mpc1b = mpcb;
0291     mpc1b.bus(:, PD) = mpc1b.bus(:, PD) .* rb.bus(:, VM).^2;
0292     mpc1b.bus(:, QD) = mpc1b.bus(:, QD) .* rb.bus(:, VM).^2;
0293     mpc1t = mpct;
0294     mpc1t.bus(:, PD) = mpc1t.bus(:, PD) .* r.bus(:, VM).^2;
0295     mpc1t.bus(:, QD) = mpc1t.bus(:, QD) .* r.bus(:, VM).^2;
0296     r1 = runcpf(mpc1b, mpc1t, mpopt0);
0297     t_ok(r1.success, [t 'success']);
0298     t_is(r1.cpf.iterations, iterations(k,3), 12, [t 'iterations']);
0299     t_ok(r.success, [t 'success']);
0300     t_is(r.bus(:, [VM VA]), r1.bus(:, [VM VA]), 6, [t 'bus voltages']);
0301     t_is(tPd, r1.bus(:, PD), 6, [t 'bus P loads']);
0302     t_is(tPd, [0; 0; 0; 0; 132.977098; 0; 154.969557; 0; 172.896177], 6, [t 'bus P loads']);
0303     t_is(tQd, r1.bus(:, QD), 6, [t 'bus Q loads']);
0304     t_is(tQd, [0; 0; 0; 0; 44.325699; 0; 54.239345; 0; 69.158471], 6, [t 'bus Q loads']);
0305     t_is(r.gen(:, [PG QG]), r1.gen(:, [PG QG]), 6, [t 'gen dispatches']);
0306     t_is(r.branch, r1.branch, 6, [t 'branch']);
0307 
0308     t = sprintf('CPF (%s) - combo ZIP loads : ', cpfparm{k});
0309     zip_w = [0.1 0.4 0.5];
0310     mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads.pw', zip_w);
0311     rb = runpf(mpcb, mpopt);
0312     r = runcpf(mpcb, mpct, mpopt);
0313     [tPd, tQd] = total_load(r, 'bus', [], mpopt);
0314     Vm = abs(rb.bus(:, VM));
0315     scaleb = [Vm.^0 Vm Vm.^2] * zip_w';
0316     Vm = abs(r.bus(:, VM));
0317     scalet = [Vm.^0 Vm Vm.^2] * zip_w';
0318     mpc1b = mpcb;
0319     mpc1b.bus(:, PD) = mpc1b.bus(:, PD) .* scaleb;
0320     mpc1b.bus(:, QD) = mpc1b.bus(:, QD) .* scaleb;
0321     mpc1t = mpct;
0322     mpc1t.bus(:, PD) = mpc1t.bus(:, PD) .* scalet;
0323     mpc1t.bus(:, QD) = mpc1t.bus(:, QD) .* scalet;
0324     r1 = runcpf(mpc1b, mpc1t, mpopt0);
0325     t_ok(r1.success, [t 'success']);
0326     t_is(r1.cpf.iterations, iterations(k,4), 12, [t 'iterations']);
0327     t_ok(r.success, [t 'success']);
0328     t_is(r.bus(:, [VM VA]), r1.bus(:, [VM VA]), 6, [t 'bus voltages']);
0329     t_is(tPd, r1.bus(:, PD), 6, [t 'bus P loads']);
0330     t_is(tPd, [0; 0; 0; 0; 139.202168; 0; 160.005294; 0; 184.209518], 6, [t 'bus P loads']);
0331     t_is(tQd, r1.bus(:, QD), 6, [t 'bus Q loads']);
0332     t_is(tQd, [0; 0; 0; 0; 46.400723; 0; 56.001853; 0; 73.683807], 6, [t 'bus Q loads']);
0333     t_is(r.gen(:, [PG QG]), r1.gen(:, [PG QG]), 6, [t 'gen dispatches']);
0334     t_is(r.branch, r1.branch, 6, [t 'branch']);
0335 end
0336 
0337 
0338 %% set up base and target cases
0339     mpcb = loadcase(casefile);
0340     mpct = mpcb;
0341     factor = 2.5 * 0.6;
0342     mpct.gen(:, [PG QG]) = mpct.gen(:, [PG QG]) * factor;
0343     mpct.bus(:, [PD QD]) = mpct.bus(:, [PD QD]) * factor;
0344     mpct.gen(:, QMAX) = 75;
0345     mpcb.gen(:, QMAX) = 75;
0346 
0347     %% with Q limits enforced
0348     iterations = [
0349         22 20 20 20;
0350         22 21 21 21;
0351         22 21 21 21;
0352     ];
0353 
0354 for k = 1:length(cpfparm)
0355     mpopt0 = mpoption(mpopt0, 'cpf.parameterization', k, 'cpf.enforce_q_lims', 1);
0356 
0357     t = sprintf('CPF (%s) w/Q lims - base case : ', cpfparm{k});
0358     r0 = runcpf(mpcb, mpct, mpopt0);
0359     t_ok(r0.success, [t 'success']);
0360     t_is(r0.cpf.iterations, iterations(k, 1), 12, [t 'iterations']);
0361     t_is(r0.cpf.max_lam, 1, 12, [t 'max_lam']);
0362 
0363     mpopt = mpopt0;
0364     t = sprintf('CPF (%s) w/Q lims - constant power only : ', cpfparm{k});
0365     zip_w = [1 0 0];
0366     mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads.pw', zip_w);
0367     r = runcpf(mpcb, mpct, mpopt);
0368     t_ok(r.success, [t 'success']);
0369     t_is(r.cpf.iterations, r0.cpf.iterations, 12, [t 'iterations']);
0370     t_is(r.cpf.max_lam, r0.cpf.max_lam, 12, [t 'max_lam']);
0371     t_is(r.bus, r0.bus, 6, [t 'bus']);
0372     t_is(r.gen, r0.gen, 6, [t 'gen']);
0373     t_is(r.branch, r0.branch, 6, [t 'branch']);
0374 
0375     t = sprintf('CPF (%s) w/Q lims - constant current only : ', cpfparm{k});
0376     zip_w = [0 1 0];
0377     mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads.pw', zip_w);
0378     rb = runpf(mpcb, mpopt);
0379     r = runcpf(mpcb, mpct, mpopt);
0380     [tPd, tQd] = total_load(r, 'bus', [], mpopt);
0381     mpc1b = mpcb;
0382     mpc1b.bus(:, PD) = mpc1b.bus(:, PD) .* rb.bus(:, VM);
0383     mpc1b.bus(:, QD) = mpc1b.bus(:, QD) .* rb.bus(:, VM);
0384     mpc1t = mpct;
0385     mpc1t.bus(:, PD) = mpc1t.bus(:, PD) .* r.bus(:, VM);
0386     mpc1t.bus(:, QD) = mpc1t.bus(:, QD) .* r.bus(:, VM);
0387     r1 = runcpf(mpc1b, mpc1t, mpopt0);
0388     t_ok(r1.success, [t 'success']);
0389     t_is(r1.cpf.iterations, iterations(k,2), 12, [t 'iterations']);
0390     t_ok(r.success, [t 'success']);
0391     t_is(r.bus(:, [VM VA]), r1.bus(:, [VM VA]), 6, [t 'bus voltages']);
0392     t_is(tPd, r1.bus(:, PD), 6, [t 'bus P loads']);
0393     t_is(tPd, [0; 0; 0; 0; 126.246737; 0; 143.114682; 0; 170.231775], 6, [t 'bus P loads']);
0394     t_is(tQd, r1.bus(:, QD), 6, [t 'bus Q loads']);
0395     t_is(tQd, [0; 0; 0; 0; 42.082245; 0; 50.090138; 0; 68.092710], 6, [t 'bus Q loads']);
0396     t_is(r.gen(:, [PG QG]), r1.gen(:, [PG QG]), 6, [t 'gen dispatches']);
0397     t_is(r.branch, r1.branch, 6, [t 'branch']);
0398 
0399     t = sprintf('CPF (%s) w/Q lims - constant impedance only : ', cpfparm{k});
0400     zip_w = [0 0 1];
0401     mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads.pw', zip_w);
0402     rb = runpf(mpcb, mpopt);
0403     r = runcpf(mpcb, mpct, mpopt);
0404     [tPd, tQd] = total_load(r, 'bus', [], mpopt);
0405     mpc1b = mpcb;
0406     mpc1b.bus(:, PD) = mpc1b.bus(:, PD) .* rb.bus(:, VM).^2;
0407     mpc1b.bus(:, QD) = mpc1b.bus(:, QD) .* rb.bus(:, VM).^2;
0408     mpc1t = mpct;
0409     mpc1t.bus(:, PD) = mpc1t.bus(:, PD) .* r.bus(:, VM).^2;
0410     mpc1t.bus(:, QD) = mpc1t.bus(:, QD) .* r.bus(:, VM).^2;
0411     r1 = runcpf(mpc1b, mpc1t, mpopt0);
0412     t_ok(r1.success, [t 'success']);
0413     t_is(r1.cpf.iterations, iterations(k,3), 12, [t 'iterations']);
0414     t_ok(r.success, [t 'success']);
0415     t_is(r.bus(:, [VM VA]), r1.bus(:, [VM VA]), 6, [t 'bus voltages']);
0416     t_is(tPd, r1.bus(:, PD), 6, [t 'bus P loads']);
0417     t_is(tPd, [0; 0; 0; 0; 119.354069; 0; 137.673893; 0; 157.193899], 6, [t 'bus P loads']);
0418     t_is(tQd, r1.bus(:, QD), 6, [t 'bus Q loads']);
0419     t_is(tQd, [0; 0; 0; 0; 39.784689; 0; 48.185862; 0; 62.877559], 6, [t 'bus Q loads']);
0420     t_is(r.gen(:, [PG QG]), r1.gen(:, [PG QG]), 6, [t 'gen dispatches']);
0421     t_is(r.branch, r1.branch, 6, [t 'branch']);
0422 
0423     t = sprintf('CPF (%s) w/Q lims - combo ZIP loads : ', cpfparm{k});
0424     zip_w = [0.1 0.4 0.5];
0425     mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads.pw', zip_w);
0426     rb = runpf(mpcb, mpopt);
0427     r = runcpf(mpcb, mpct, mpopt);
0428     [tPd, tQd] = total_load(r, 'bus', [], mpopt);
0429     Vm = abs(rb.bus(:, VM));
0430     scaleb = [Vm.^0 Vm Vm.^2] * zip_w';
0431     Vm = abs(r.bus(:, VM));
0432     scalet = [Vm.^0 Vm Vm.^2] * zip_w';
0433     mpc1b = mpcb;
0434     mpc1b.bus(:, PD) = mpc1b.bus(:, PD) .* scaleb;
0435     mpc1b.bus(:, QD) = mpc1b.bus(:, QD) .* scaleb;
0436     mpc1t = mpct;
0437     mpc1t.bus(:, PD) = mpc1t.bus(:, PD) .* scalet;
0438     mpc1t.bus(:, QD) = mpc1t.bus(:, QD) .* scalet;
0439     r1 = runcpf(mpc1b, mpc1t, mpopt0);
0440     t_ok(r1.success, [t 'success']);
0441     t_is(r1.cpf.iterations, iterations(k,4), 12, [t 'iterations']);
0442     t_ok(r.success, [t 'success']);
0443     t_is(r.bus(:, [VM VA]), r1.bus(:, [VM VA]), 6, [t 'bus voltages']);
0444     t_is(tPd, r1.bus(:, PD), 6, [t 'bus P loads']);
0445     t_is(tPd, [0; 0; 0; 0; 123.418124; 0; 140.853612; 0; 164.913503], 6, [t 'bus P loads']);
0446     t_is(tQd, r1.bus(:, QD), 6, [t 'bus Q loads']);
0447     t_is(tQd, [0; 0; 0; 0; 41.1393746; 0; 49.298764; 0; 65.9654012], 6, [t 'bus Q loads']);
0448     t_is(r.gen(:, [PG QG]), r1.gen(:, [PG QG]), 6, [t 'gen dispatches']);
0449     t_is(r.branch, r1.branch, 6, [t 'branch']);
0450 end
0451 
0452 t_end;