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
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0067 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0068 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0069 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0070 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0071 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0072 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0073 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0074 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0075
0076
0077 alg = upper(mpopt.opf.dc.solver);
0078
0079
0080 mpc = get_mpc(om);
0081 [baseMVA, bus, gen, branch, gencost] = ...
0082 deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost);
0083 cp = get_cost_params(om);
0084 [N, H, Cw] = deal(cp.N, cp.H, cp.Cw);
0085 fparm = [cp.dd cp.rh cp.kk cp.mm];
0086 Bf = userdata(om, 'Bf');
0087 Pfinj = userdata(om, 'Pfinj');
0088 [vv, ll] = get_idx(om);
0089
0090
0091 ipol = find(gencost(:, MODEL) == POLYNOMIAL);
0092 ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0093 nb = size(bus, 1);
0094 nl = size(branch, 1);
0095 nw = size(N, 1);
0096 ny = getN(om, 'var', 'y');
0097 nxyz = getN(om, 'var');
0098
0099
0100 [A, l, u] = linear_constraints(om);
0101 [x0, xmin, xmax] = getv(om);
0102
0103
0104
0105
0106
0107
0108
0109 any_pwl = (ny > 0);
0110 if any_pwl
0111 Npwl = sparse(ones(ny,1), vv.i1.y:vv.iN.y, 1, 1, nxyz);
0112 Hpwl = 0;
0113 Cpwl = 1;
0114 fparm_pwl = [1 0 0 1];
0115 else
0116 Npwl = sparse(0, nxyz);
0117 Hpwl = [];
0118 Cpwl = [];
0119 fparm_pwl = [];
0120 end
0121
0122
0123 npol = length(ipol);
0124 if any(find(gencost(ipol, NCOST) > 3))
0125 error('DC opf cannot handle polynomial costs with higher than quadratic order.');
0126 end
0127 iqdr = find(gencost(ipol, NCOST) == 3);
0128 ilin = find(gencost(ipol, NCOST) == 2);
0129 polycf = zeros(npol, 3);
0130 if ~isempty(iqdr)
0131 polycf(iqdr, :) = gencost(ipol(iqdr), COST:COST+2);
0132 end
0133 polycf(ilin, 2:3) = gencost(ipol(ilin), COST:COST+1);
0134 polycf = polycf * diag([ baseMVA^2 baseMVA 1]);
0135 Npol = sparse(1:npol, vv.i1.Pg-1+ipol, 1, npol, nxyz);
0136 Hpol = sparse(1:npol, 1:npol, 2*polycf(:, 1), npol, npol);
0137 Cpol = polycf(:, 2);
0138 fparm_pol = ones(npol,1) * [ 1 0 0 1 ];
0139
0140
0141 NN = [ Npwl; Npol; N ];
0142 HHw = [ Hpwl, sparse(any_pwl, npol+nw);
0143 sparse(npol, any_pwl), Hpol, sparse(npol, nw);
0144 sparse(nw, any_pwl+npol), H ];
0145 CCw = [Cpwl; Cpol; Cw];
0146 ffparm = [ fparm_pwl; fparm_pol; fparm ];
0147
0148
0149 nnw = any_pwl+npol+nw;
0150 M = sparse(1:nnw, 1:nnw, ffparm(:, 4), nnw, nnw);
0151 MR = M * ffparm(:, 2);
0152 HMR = HHw * MR;
0153 MN = M * NN;
0154 HH = MN' * HHw * MN;
0155 CC = full(MN' * (CCw - HMR));
0156 C0 = 1/2 * MR' * HMR + sum(polycf(:, 3));
0157
0158
0159 if strcmp(alg, 'DEFAULT')
0160 if have_fcn('cplex')
0161 alg = 'CPLEX';
0162 elseif have_fcn('gurobi')
0163 alg = 'GUROBI';
0164 elseif have_fcn('mosek')
0165 alg = 'MOSEK';
0166 elseif have_fcn('bpmpd')
0167 alg = 'BPMPD';
0168 elseif have_fcn('quadprog')
0169 alg = 'OT';
0170 elseif (isempty(HH) || ~any(any(HH))) && have_fcn('glpk')
0171 alg = 'GLPK';
0172 else
0173 alg = 'MIPS';
0174 end
0175 end
0176
0177
0178 opt = struct('alg', alg, 'verbose', mpopt.verbose);
0179 switch alg
0180 case {'MIPS', 'IPOPT'}
0181
0182 if mpopt.opf.init_from_mpc ~= 1
0183 Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180);
0184
0185 lb = xmin; ub = xmax;
0186 lb(xmin == -Inf) = -1e10;
0187 ub(xmax == Inf) = 1e10;
0188 x0 = (lb + ub) / 2;
0189 k = find(xmin == -Inf & xmax < Inf);
0190 x0(k) = xmax(k) - 1;
0191 k = find(xmin > -Inf & xmax == Inf);
0192 x0(k) = xmin(k) + 1;
0193 x0(vv.i1.Va:vv.iN.Va) = Varefs(1);
0194 if ny > 0
0195 ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0196 c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST)));
0197 x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c));
0198 end
0199 end
0200 switch alg
0201 case 'MIPS'
0202
0203 opt.mips_opt = mpopt.mips;
0204 opt.mips_opt.verbose = mpopt.verbose;
0205 if opt.mips_opt.feastol == 0
0206 opt.mips_opt.feastol = mpopt.opf.violation;
0207 end
0208 if ~isfield(opt.mips_opt, 'cost_mult') || isempty(opt.mips_opt.cost_mult)
0209 opt.mips_opt.cost_mult = 1;
0210 end
0211 case 'IPOPT'
0212 opt.ipopt_opt = ipopt_options([], mpopt);
0213 end
0214 case 'CPLEX'
0215 opt.cplex_opt = cplex_options([], mpopt);
0216 case 'GLPK'
0217 opt.glpk_opt = glpk_options([], mpopt);
0218 case 'GUROBI'
0219 opt.grb_opt = gurobi_options([], mpopt);
0220 case 'MOSEK'
0221 opt.mosek_opt = mosek_options([], mpopt);
0222 case 'OT'
0223 if isfield(mpopt, 'linprog') && ~isempty(mpopt.linprog)
0224 opt.linprog_opt = mpopt.linprog;
0225 end
0226 if isfield(mpopt, 'quadprog') && ~isempty(mpopt.quadprog)
0227 opt.quadprog_opt = mpopt.quadprog;
0228 end
0229 end
0230
0231
0232 [x, f, info, output, lambda] = qps_matpower(HH, CC, A, l, u, xmin, xmax, x0, opt);
0233 success = (info == 1);
0234
0235
0236 if ~any(isnan(x))
0237
0238 Va = x(vv.i1.Va:vv.iN.Va);
0239 Pg = x(vv.i1.Pg:vv.iN.Pg);
0240 f = f + C0;
0241
0242
0243 bus(:, VM) = ones(nb, 1);
0244 bus(:, VA) = Va * 180/pi;
0245 gen(:, PG) = Pg * baseMVA;
0246
0247
0248 branch(:, [QF, QT]) = zeros(nl, 2);
0249 branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;
0250 branch(:, PT) = -branch(:, PF);
0251 end
0252
0253
0254 mu_l = lambda.mu_l;
0255 mu_u = lambda.mu_u;
0256 muLB = lambda.lower;
0257 muUB = lambda.upper;
0258
0259
0260 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0261 bus(:, [LAM_P, LAM_Q, MU_VMIN, MU_VMAX]) = zeros(nb, 4);
0262 gen(:, [MU_PMIN, MU_PMAX, MU_QMIN, MU_QMAX]) = zeros(size(gen, 1), 4);
0263 branch(:, [MU_SF, MU_ST]) = zeros(nl, 2);
0264 bus(:, LAM_P) = (mu_u(ll.i1.Pmis:ll.iN.Pmis) - mu_l(ll.i1.Pmis:ll.iN.Pmis)) / baseMVA;
0265 branch(il, MU_SF) = mu_u(ll.i1.Pf:ll.iN.Pf) / baseMVA;
0266 branch(il, MU_ST) = mu_l(ll.i1.Pf:ll.iN.Pf) / baseMVA;
0267 gen(:, MU_PMIN) = muLB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0268 gen(:, MU_PMAX) = muUB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0269 pimul = [
0270 mu_l - mu_u;
0271 -ones(ny>0, 1);
0272 muLB - muUB
0273 ];
0274
0275 mu = struct( ...
0276 'var', struct('l', muLB, 'u', muUB), ...
0277 'lin', struct('l', mu_l, 'u', mu_u) );
0278
0279 results = mpc;
0280 [results.bus, results.branch, results.gen, ...
0281 results.om, results.x, results.mu, results.f] = ...
0282 deal(bus, branch, gen, om, x, mu, f);
0283
0284 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', output);