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