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