0001 function [results, success, raw] = mipsopf_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
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0070 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0071 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0072 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0073 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0074 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0075 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0076 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0077 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0078
0079
0080 opt = mpopt.mips;
0081 opt.verbose = mpopt.verbose;
0082 if opt.feastol == 0
0083 opt.feastol = mpopt.opf.violation;
0084 end
0085 if ~isfield(opt, 'cost_mult') || isempty(opt.cost_mult)
0086 opt.cost_mult = 1e-4;
0087 end
0088
0089
0090 mpc = get_mpc(om);
0091 [baseMVA, bus, gen, branch, gencost] = ...
0092 deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0093 [vv, ll, nn] = get_idx(om);
0094
0095
0096 nb = size(bus, 1);
0097 nl = size(branch, 1);
0098 ny = getN(om, 'var', 'y');
0099
0100
0101 [A, l, u] = linear_constraints(om);
0102
0103
0104 [x0, xmin, xmax] = getv(om);
0105
0106
0107 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0108
0109
0110 if mpopt.opf.init_from_mpc ~= 1
0111 ll = xmin; uu = xmax;
0112 ll(xmin == -Inf) = -1e10;
0113 uu(xmax == Inf) = 1e10;
0114 x0 = (ll + uu) / 2;
0115 k = find(xmin == -Inf & xmax < Inf);
0116 x0(k) = xmax(k) - 1;
0117 k = find(xmin > -Inf & xmax == Inf);
0118 x0(k) = xmin(k) + 1;
0119 Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0120 x0(vv.i1.Va:vv.iN.Va) = Varefs(1);
0121 if ny > 0
0122 ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0123
0124
0125 c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));
0126 x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0127
0128 end
0129 end
0130
0131
0132 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0133 nl2 = length(il);
0134
0135
0136 f_fcn = @(x)opf_costfcn(x, om);
0137 gh_fcn = @(x)opf_consfcn(x, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0138 hess_fcn = @(x, lambda, cost_mult)opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0139 [x, f, info, Output, Lambda] = ...
0140 mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);
0141 success = (info > 0);
0142
0143
0144 Va = x(vv.i1.Va:vv.iN.Va);
0145 Vm = x(vv.i1.Vm:vv.iN.Vm);
0146 Pg = x(vv.i1.Pg:vv.iN.Pg);
0147 Qg = x(vv.i1.Qg:vv.iN.Qg);
0148 V = Vm .* exp(1j*Va);
0149
0150
0151
0152 bus(:, VA) = Va * 180/pi;
0153 bus(:, VM) = Vm;
0154 gen(:, PG) = Pg * baseMVA;
0155 gen(:, QG) = Qg * baseMVA;
0156 gen(:, VG) = Vm(gen(:, GEN_BUS));
0157
0158
0159 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);
0160 St = V(branch(:, T_BUS)) .* conj(Yt * V);
0161 branch(:, PF) = real(Sf) * baseMVA;
0162 branch(:, QF) = imag(Sf) * baseMVA;
0163 branch(:, PT) = real(St) * baseMVA;
0164 branch(:, QT) = imag(St) * baseMVA;
0165
0166
0167
0168 muSf = zeros(nl, 1);
0169 muSt = zeros(nl, 1);
0170 if ~isempty(il)
0171 muSf(il) = 2 * Lambda.ineqnonlin(1:nl2) .* branch(il, RATE_A) / baseMVA;
0172 muSt(il) = 2 * Lambda.ineqnonlin((1:nl2)+nl2) .* branch(il, RATE_A) / baseMVA;
0173 end
0174
0175
0176 bus(:, MU_VMAX) = Lambda.upper(vv.i1.Vm:vv.iN.Vm);
0177 bus(:, MU_VMIN) = Lambda.lower(vv.i1.Vm:vv.iN.Vm);
0178 gen(:, MU_PMAX) = Lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0179 gen(:, MU_PMIN) = Lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0180 gen(:, MU_QMAX) = Lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0181 gen(:, MU_QMIN) = Lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0182 bus(:, LAM_P) = Lambda.eqnonlin(nn.i1.Pmis:nn.iN.Pmis) / baseMVA;
0183 bus(:, LAM_Q) = Lambda.eqnonlin(nn.i1.Qmis:nn.iN.Qmis) / baseMVA;
0184 branch(:, MU_SF) = muSf / baseMVA;
0185 branch(:, MU_ST) = muSt / baseMVA;
0186
0187
0188 nlnN = getN(om, 'nln');
0189
0190
0191 kl = find(Lambda.eqnonlin < 0);
0192 ku = find(Lambda.eqnonlin > 0);
0193 nl_mu_l = zeros(nlnN, 1);
0194 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0195 nl_mu_l(kl) = -Lambda.eqnonlin(kl);
0196 nl_mu_u(ku) = Lambda.eqnonlin(ku);
0197
0198 mu = struct( ...
0199 'var', struct('l', Lambda.lower, 'u', Lambda.upper), ...
0200 'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0201 'lin', struct('l', Lambda.mu_l, 'u', Lambda.mu_u) );
0202
0203 results = mpc;
0204 [results.bus, results.branch, results.gen, ...
0205 results.om, results.x, results.mu, results.f] = ...
0206 deal(bus, branch, gen, om, x, mu, f);
0207
0208 pimul = [ ...
0209 results.mu.nln.l - results.mu.nln.u;
0210 results.mu.lin.l - results.mu.lin.u;
0211 -ones(ny>0, 1);
0212 results.mu.var.l - results.mu.var.u;
0213 ];
0214 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', Output);