0001 function [results, success, raw] = nlpopf_solver(om, mpopt)
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
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0051 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0052 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0053 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0054 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0055 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0056 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0057 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0058 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072 mpc = om.get_mpc();
0073 [baseMVA, bus, gen, branch, gencost] = ...
0074 deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0075 [vv, ll, nne, nni] = om.get_idx();
0076
0077
0078 nb = size(bus, 1);
0079 nl = size(branch, 1);
0080 ny = om.getN('var', 'y');
0081
0082
0083 model = om.problem_type();
0084 opt = mpopt2nlpopt(mpopt, model);
0085
0086
0087 if mpopt.opf.start < 2
0088 [x0, xmin, xmax] = om.params_var();
0089 s = 1;
0090 lb = xmin; ub = xmax;
0091 lb(xmin == -Inf) = -1e10;
0092 ub(xmax == Inf) = 1e10;
0093 x0 = (lb + ub) / 2;
0094 k = find(xmin == -Inf & xmax < Inf);
0095 x0(k) = xmax(k) - s;
0096 k = find(xmin > -Inf & xmax == Inf);
0097 x0(k) = xmin(k) + s;
0098 Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0099 Vmax = min(bus(:, VMAX), 1.5);
0100 Vmin = max(bus(:, VMIN), 0.5);
0101 Vm = (Vmax + Vmin) / 2;
0102 if mpopt.opf.v_cartesian
0103 V = Vm * exp(1j*Varefs(1));
0104 x0(vv.i1.Vr:vv.iN.Vr) = real(V);
0105 x0(vv.i1.Vi:vv.iN.Vi) = imag(V);
0106 else
0107 x0(vv.i1.Va:vv.iN.Va) = Varefs(1);
0108 x0(vv.i1.Vm:vv.iN.Vm) = Vm;
0109 if ny > 0
0110 ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0111 c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));
0112 x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0113 end
0114 end
0115 opt.x0 = x0;
0116 end
0117
0118
0119 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0120 nl2 = length(il);
0121
0122
0123 [x, f, eflag, output, lambda] = om.solve(opt);
0124 success = (eflag > 0);
0125
0126
0127 if mpopt.opf.v_cartesian
0128 Vi = x(vv.i1.Vi:vv.iN.Vi);
0129 Vr = x(vv.i1.Vr:vv.iN.Vr);
0130 V = Vr + 1j*Vi;
0131 Va = angle(V);
0132 Vm = abs(V);
0133 else
0134 Va = x(vv.i1.Va:vv.iN.Va);
0135 Vm = x(vv.i1.Vm:vv.iN.Vm);
0136 V = Vm .* exp(1j*Va);
0137 end
0138 Pg = x(vv.i1.Pg:vv.iN.Pg);
0139 Qg = x(vv.i1.Qg:vv.iN.Qg);
0140
0141
0142
0143 bus(:, VA) = Va * 180/pi;
0144 bus(:, VM) = Vm;
0145 gen(:, PG) = Pg * baseMVA;
0146 gen(:, QG) = Qg * baseMVA;
0147 gen(:, VG) = Vm(gen(:, GEN_BUS));
0148
0149
0150 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0151 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);
0152 St = V(branch(:, T_BUS)) .* conj(Yt * V);
0153 branch(:, PF) = real(Sf) * baseMVA;
0154 branch(:, QF) = imag(Sf) * baseMVA;
0155 branch(:, PT) = real(St) * baseMVA;
0156 branch(:, QT) = imag(St) * baseMVA;
0157
0158
0159
0160 muSf = zeros(nl, 1);
0161 muSt = zeros(nl, 1);
0162 if ~isempty(il)
0163 if upper(mpopt.opf.flow_lim(1)) == 'P'
0164 muSf(il) = lambda.ineqnonlin(nni.i1.Sf:nni.iN.Sf);
0165 muSt(il) = lambda.ineqnonlin(nni.i1.St:nni.iN.St);
0166 else
0167 muSf(il) = 2 * lambda.ineqnonlin(nni.i1.Sf:nni.iN.Sf) .* branch(il, RATE_A) / baseMVA;
0168 muSt(il) = 2 * lambda.ineqnonlin(nni.i1.St:nni.iN.St) .* branch(il, RATE_A) / baseMVA;
0169 end
0170 end
0171
0172
0173 if mpopt.opf.v_cartesian
0174 if om.userdata.veq
0175 lam = lambda.eqnonlin(nne.i1.Veq:nne.iN.Veq);
0176 mu_Vmax = zeros(size(lam));
0177 mu_Vmin = zeros(size(lam));
0178 mu_Vmax(lam > 0) = lam(lam > 0);
0179 mu_Vmin(lam < 0) = -lam(lam < 0);
0180 bus(om.userdata.veq, MU_VMAX) = mu_Vmax;
0181 bus(om.userdata.veq, MU_VMIN) = mu_Vmin;
0182 end
0183 bus(om.userdata.viq, MU_VMAX) = lambda.ineqnonlin(nni.i1.Vmax:nni.iN.Vmax);
0184 bus(om.userdata.viq, MU_VMIN) = lambda.ineqnonlin(nni.i1.Vmin:nni.iN.Vmin);
0185 else
0186 bus(:, MU_VMAX) = lambda.upper(vv.i1.Vm:vv.iN.Vm);
0187 bus(:, MU_VMIN) = lambda.lower(vv.i1.Vm:vv.iN.Vm);
0188 end
0189 gen(:, MU_PMAX) = lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0190 gen(:, MU_PMIN) = lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0191 gen(:, MU_QMAX) = lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0192 gen(:, MU_QMIN) = lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0193 if mpopt.opf.current_balance
0194
0195
0196
0197
0198
0199
0200 VV = V ./ (V .* conj(V));
0201 VVr = real(VV);
0202 VVi = imag(VV);
0203 lamM = lambda.eqnonlin(nne.i1.rImis:nne.iN.rImis);
0204 lamN = lambda.eqnonlin(nne.i1.iImis:nne.iN.iImis);
0205 bus(:, LAM_P) = (VVr.*lamM + VVi.*lamN) / baseMVA;
0206 bus(:, LAM_Q) = (VVi.*lamM - VVr.*lamN) / baseMVA;
0207 else
0208 bus(:, LAM_P) = lambda.eqnonlin(nne.i1.Pmis:nne.iN.Pmis) / baseMVA;
0209 bus(:, LAM_Q) = lambda.eqnonlin(nne.i1.Qmis:nne.iN.Qmis) / baseMVA;
0210 end
0211 branch(:, MU_SF) = muSf / baseMVA;
0212 branch(:, MU_ST) = muSt / baseMVA;
0213
0214
0215 nlnN = 2*nb + 2*nl;
0216
0217
0218 kl = find(lambda.eqnonlin(1:2*nb) < 0);
0219 ku = find(lambda.eqnonlin(1:2*nb) > 0);
0220 nl_mu_l = zeros(nlnN, 1);
0221 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0222 nl_mu_l(kl) = -lambda.eqnonlin(kl);
0223 nl_mu_u(ku) = lambda.eqnonlin(ku);
0224
0225 if isfield(lambda, 'ineqnonlin')
0226 lam_nli = lambda.ineqnonlin;
0227 else
0228 lam_nli = [];
0229 end
0230
0231 mu = struct( ...
0232 'var', struct('l', lambda.lower, 'u', lambda.upper), ...
0233 'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0234 'nle', lambda.eqnonlin, ...
0235 'nli', lam_nli, ...
0236 'lin', struct('l', lambda.mu_l, 'u', lambda.mu_u) );
0237
0238 results = mpc;
0239 [results.bus, results.branch, results.gen, ...
0240 results.om, results.x, results.mu, results.f] = ...
0241 deal(bus, branch, gen, om, x, mu, f);
0242
0243 pimul = [ ...
0244 results.mu.nln.l - results.mu.nln.u;
0245 results.mu.lin.l - results.mu.lin.u;
0246 -ones(ny>0, 1);
0247 results.mu.var.l - results.mu.var.u;
0248 ];
0249 raw = struct('xr', x, 'pimul', pimul, 'info', eflag, 'output', output);