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 [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 opt = mpopt.mips;
0062 opt.verbose = mpopt.verbose;
0063 if opt.feastol == 0
0064 opt.feastol = mpopt.opf.violation;
0065 end
0066 if ~isfield(opt, 'cost_mult') || isempty(opt.cost_mult)
0067 opt.cost_mult = 1e-4;
0068 end
0069
0070
0071 mpc = get_mpc(om);
0072 [baseMVA, bus, gen, branch, gencost] = ...
0073 deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0074 [vv, ll, nn] = get_idx(om);
0075
0076
0077 nb = size(bus, 1);
0078 nl = size(branch, 1);
0079 ny = getN(om, 'var', 'y');
0080
0081
0082 [A, l, u] = linear_constraints(om);
0083
0084
0085 [x0, xmin, xmax] = getv(om);
0086
0087
0088 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0089
0090
0091 if mpopt.opf.init_from_mpc ~= 1
0092 ll = xmin; uu = xmax;
0093 ll(xmin == -Inf) = -1e10;
0094 uu(xmax == Inf) = 1e10;
0095 x0 = (ll + uu) / 2;
0096 k = find(xmin == -Inf & xmax < Inf);
0097 x0(k) = xmax(k) - 1;
0098 k = find(xmin > -Inf & xmax == Inf);
0099 x0(k) = xmin(k) + 1;
0100 Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0101 x0(vv.i1.Va:vv.iN.Va) = Varefs(1);
0102 if ny > 0
0103 ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0104
0105
0106 c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));
0107 x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0108
0109 end
0110 end
0111
0112
0113 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0114 nl2 = length(il);
0115
0116
0117 f_fcn = @(x)opf_costfcn(x, om);
0118 gh_fcn = @(x)opf_consfcn(x, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0119 hess_fcn = @(x, lambda, cost_mult)opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0120 [x, f, info, Output, Lambda] = ...
0121 mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);
0122 success = (info > 0);
0123
0124
0125 Va = x(vv.i1.Va:vv.iN.Va);
0126 Vm = x(vv.i1.Vm:vv.iN.Vm);
0127 Pg = x(vv.i1.Pg:vv.iN.Pg);
0128 Qg = x(vv.i1.Qg:vv.iN.Qg);
0129 V = Vm .* exp(1j*Va);
0130
0131
0132
0133 bus(:, VA) = Va * 180/pi;
0134 bus(:, VM) = Vm;
0135 gen(:, PG) = Pg * baseMVA;
0136 gen(:, QG) = Qg * baseMVA;
0137 gen(:, VG) = Vm(gen(:, GEN_BUS));
0138
0139
0140 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);
0141 St = V(branch(:, T_BUS)) .* conj(Yt * V);
0142 branch(:, PF) = real(Sf) * baseMVA;
0143 branch(:, QF) = imag(Sf) * baseMVA;
0144 branch(:, PT) = real(St) * baseMVA;
0145 branch(:, QT) = imag(St) * baseMVA;
0146
0147
0148
0149 muSf = zeros(nl, 1);
0150 muSt = zeros(nl, 1);
0151 if ~isempty(il)
0152 muSf(il) = 2 * Lambda.ineqnonlin(1:nl2) .* branch(il, RATE_A) / baseMVA;
0153 muSt(il) = 2 * Lambda.ineqnonlin((1:nl2)+nl2) .* branch(il, RATE_A) / baseMVA;
0154 end
0155
0156
0157 bus(:, MU_VMAX) = Lambda.upper(vv.i1.Vm:vv.iN.Vm);
0158 bus(:, MU_VMIN) = Lambda.lower(vv.i1.Vm:vv.iN.Vm);
0159 gen(:, MU_PMAX) = Lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0160 gen(:, MU_PMIN) = Lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0161 gen(:, MU_QMAX) = Lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0162 gen(:, MU_QMIN) = Lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0163 bus(:, LAM_P) = Lambda.eqnonlin(nn.i1.Pmis:nn.iN.Pmis) / baseMVA;
0164 bus(:, LAM_Q) = Lambda.eqnonlin(nn.i1.Qmis:nn.iN.Qmis) / baseMVA;
0165 branch(:, MU_SF) = muSf / baseMVA;
0166 branch(:, MU_ST) = muSt / baseMVA;
0167
0168
0169 nlnN = getN(om, 'nln');
0170
0171
0172 kl = find(Lambda.eqnonlin < 0);
0173 ku = find(Lambda.eqnonlin > 0);
0174 nl_mu_l = zeros(nlnN, 1);
0175 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0176 nl_mu_l(kl) = -Lambda.eqnonlin(kl);
0177 nl_mu_u(ku) = Lambda.eqnonlin(ku);
0178
0179 mu = struct( ...
0180 'var', struct('l', Lambda.lower, 'u', Lambda.upper), ...
0181 'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0182 'lin', struct('l', Lambda.mu_l, 'u', Lambda.mu_u) );
0183
0184 results = mpc;
0185 [results.bus, results.branch, results.gen, ...
0186 results.om, results.x, results.mu, results.f] = ...
0187 deal(bus, branch, gen, om, x, mu, f);
0188
0189 pimul = [ ...
0190 results.mu.nln.l - results.mu.nln.u;
0191 results.mu.lin.l - results.mu.lin.u;
0192 -ones(ny>0, 1);
0193 results.mu.var.l - results.mu.var.u;
0194 ];
0195 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', Output);