0001 function [results, success, raw] = dcopf_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 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0046 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0047 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0048 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0049 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0050 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0051 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0052 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0053 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0054
0055
0056 mpc = get_mpc(om);
0057 [baseMVA, bus, gen, branch, gencost] = ...
0058 deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0059 cp = get_cost_params(om);
0060 [N, H, Cw] = deal(cp.N, cp.H, cp.Cw);
0061 fparm = [cp.dd cp.rh cp.kk cp.mm];
0062 Bf = userdata(om, 'Bf');
0063 Pfinj = userdata(om, 'Pfinj');
0064 [vv, ll] = get_idx(om);
0065
0066
0067 ipol = find(gencost(:, MODEL) == POLYNOMIAL);
0068 ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0069 nb = size(bus, 1);
0070 nl = size(branch, 1);
0071 nw = size(N, 1);
0072 ny = getN(om, 'var', 'y');
0073 nxyz = getN(om, 'var');
0074
0075
0076 [A, l, u] = linear_constraints(om);
0077 [x0, xmin, xmax] = getv(om);
0078
0079
0080
0081
0082
0083
0084
0085 any_pwl = (ny > 0);
0086 if any_pwl
0087 Npwl = sparse(ones(ny,1), vv.i1.y:vv.iN.y, 1, 1, nxyz);
0088 Hpwl = 0;
0089 Cpwl = 1;
0090 fparm_pwl = [1 0 0 1];
0091 else
0092 Npwl = sparse(0, nxyz);
0093 Hpwl = [];
0094 Cpwl = [];
0095 fparm_pwl = [];
0096 end
0097
0098
0099 npol = length(ipol);
0100 if any(find(gencost(ipol, NCOST) > 3))
0101 error('DC opf cannot handle polynomial costs with higher than quadratic order.');
0102 end
0103 iqdr = find(gencost(ipol, NCOST) == 3);
0104 ilin = find(gencost(ipol, NCOST) == 2);
0105 polycf = zeros(npol, 3);
0106 if ~isempty(iqdr)
0107 polycf(iqdr, :) = gencost(ipol(iqdr), COST:COST+2);
0108 end
0109 polycf(ilin, 2:3) = gencost(ipol(ilin), COST:COST+1);
0110 polycf = polycf * diag([ baseMVA^2 baseMVA 1]);
0111 Npol = sparse(1:npol, vv.i1.Pg-1+ipol, 1, npol, nxyz);
0112 Hpol = sparse(1:npol, 1:npol, 2*polycf(:, 1), npol, npol);
0113 Cpol = polycf(:, 2);
0114 fparm_pol = ones(npol,1) * [ 1 0 0 1 ];
0115
0116
0117 NN = [ Npwl; Npol; N ];
0118 HHw = [ Hpwl, sparse(any_pwl, npol+nw);
0119 sparse(npol, any_pwl), Hpol, sparse(npol, nw);
0120 sparse(nw, any_pwl+npol), H ];
0121 CCw = [Cpwl; Cpol; Cw];
0122 ffparm = [ fparm_pwl; fparm_pol; fparm ];
0123
0124
0125 nnw = any_pwl+npol+nw;
0126 M = sparse(1:nnw, 1:nnw, ffparm(:, 4), nnw, nnw);
0127 MR = M * ffparm(:, 2);
0128 HMR = HHw * MR;
0129 MN = M * NN;
0130 HH = MN' * HHw * MN;
0131 CC = full(MN' * (CCw - HMR));
0132 C0 = 1/2 * MR' * HMR + sum(polycf(:, 3));
0133
0134
0135 if isempty(HH) || ~any(any(HH))
0136 model = 'LP';
0137 else
0138 model = 'QP';
0139 end
0140 opt = mpopt2qpopt(mpopt, model);
0141
0142
0143 if mpopt.opf.init_from_mpc ~= 1 && ...
0144 (strcmp(opt.alg, 'MIPS') || strcmp(opt.alg, 'IPOPT'))
0145 Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0146
0147 lb = xmin; ub = xmax;
0148 lb(xmin == -Inf) = -1e10;
0149 ub(xmax == Inf) = 1e10;
0150 x0 = (lb + ub) / 2;
0151 k = find(xmin == -Inf & xmax < Inf);
0152 x0(k) = xmax(k) - 1;
0153 k = find(xmin > -Inf & xmax == Inf);
0154 x0(k) = xmin(k) + 1;
0155 x0(vv.i1.Va:vv.iN.Va) = Varefs(1);
0156 if ny > 0
0157 ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0158 c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));
0159 x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0160 end
0161 end
0162
0163
0164 [x, f, info, output, lambda] = qps_matpower(HH, CC, A, l, u, xmin, xmax, x0, opt);
0165 success = (info == 1);
0166
0167
0168 if ~any(isnan(x))
0169
0170 Va = x(vv.i1.Va:vv.iN.Va);
0171 Pg = x(vv.i1.Pg:vv.iN.Pg);
0172 f = f + C0;
0173
0174
0175 bus(:, VM) = ones(nb, 1);
0176 bus(:, VA) = Va * 180/pi;
0177 gen(:, PG) = Pg * baseMVA;
0178
0179
0180 branch(:, [QF, QT]) = zeros(nl, 2);
0181 branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;
0182 branch(:, PT) = -branch(:, PF);
0183 end
0184
0185
0186 mu_l = lambda.mu_l;
0187 mu_u = lambda.mu_u;
0188 muLB = lambda.lower;
0189 muUB = lambda.upper;
0190
0191
0192 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0193 bus(:, [LAM_P, LAM_Q, MU_VMIN, MU_VMAX]) = zeros(nb, 4);
0194 gen(:, [MU_PMIN, MU_PMAX, MU_QMIN, MU_QMAX]) = zeros(size(gen, 1), 4);
0195 branch(:, [MU_SF, MU_ST]) = zeros(nl, 2);
0196 bus(:, LAM_P) = (mu_u(ll.i1.Pmis:ll.iN.Pmis) - mu_l(ll.i1.Pmis:ll.iN.Pmis)) / baseMVA;
0197 branch(il, MU_SF) = mu_u(ll.i1.Pf:ll.iN.Pf) / baseMVA;
0198 branch(il, MU_ST) = mu_l(ll.i1.Pf:ll.iN.Pf) / baseMVA;
0199 gen(:, MU_PMIN) = muLB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0200 gen(:, MU_PMAX) = muUB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0201 pimul = [
0202 mu_l - mu_u;
0203 -ones(ny>0, 1);
0204 muLB - muUB
0205 ];
0206
0207 mu = struct( ...
0208 'var', struct('l', muLB, 'u', muUB), ...
0209 'lin', struct('l', mu_l, 'u', mu_u) );
0210
0211 results = mpc;
0212 [results.bus, results.branch, results.gen, ...
0213 results.om, results.x, results.mu, results.f] = ...
0214 deal(bus, branch, gen, om, x, mu, f);
0215
0216 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', output);