0001 function [results, success, raw] = fmincopf_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
0067
0068
0069 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0070 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0071 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0072 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0073 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0074 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0075 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0076 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0077
0078
0079 mpc = get_mpc(om);
0080 [baseMVA, bus, gen, branch] = ...
0081 deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch);
0082 [vv, ll, nn] = get_idx(om);
0083
0084
0085 nb = size(bus, 1);
0086 nl = size(branch, 1);
0087 ny = getN(om, 'var', 'y');
0088
0089
0090 [x0, LB, UB] = getv(om);
0091
0092
0093 [A, l, u] = linear_constraints(om);
0094
0095
0096
0097 ieq = find( abs(u-l) <= eps );
0098 igt = find( u >= 1e10 & l > -1e10 );
0099 ilt = find( l <= -1e10 & u < 1e10 );
0100 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0101 Af = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0102 bf = [ u(ilt); -l(igt); u(ibx); -l(ibx)];
0103 Afeq = A(ieq, :);
0104 bfeq = u(ieq);
0105
0106
0107 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0108
0109
0110 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0111 nl2 = length(il);
0112
0113
0114 fmoptions = optimset('GradObj', 'on', 'GradConstr', 'on', ...
0115 'TolCon', mpopt.opf.violation, 'TolX', mpopt.fmincon.tol_x, ...
0116 'TolFun', mpopt.fmincon.tol_f );
0117 if mpopt.fmincon.max_it ~= 0
0118 fmoptions = optimset(fmoptions, 'MaxIter', mpopt.fmincon.max_it, ...
0119 'MaxFunEvals', 4 * mpopt.fmincon.max_it);
0120 end
0121
0122 if mpopt.verbose == 0,
0123 fmoptions.Display = 'off';
0124 elseif mpopt.verbose == 1
0125 fmoptions.Display = 'iter';
0126 else
0127 fmoptions.Display = 'testing';
0128 end
0129
0130
0131 if have_fcn('fmincon_ipm')
0132 switch mpopt.fmincon.alg
0133 case 1
0134 fmoptions = optimset(fmoptions, 'Algorithm', 'active-set');
0135 Af = full(Af);
0136 Afeq = full(Afeq);
0137 case 2
0138 fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point');
0139 case 3
0140 fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point', 'Hessian','lbfgs');
0141 case 4
0142 fmc_hessian = @(x, lambda)opf_hessfcn(x, lambda, 1, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0143 fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point', ...
0144 'Hessian', 'user-supplied', 'HessFcn', fmc_hessian);
0145 case 5
0146 fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point', 'Hessian','fin-diff-grads', 'SubProblem', 'cg');
0147 case 6
0148 fmoptions = optimset(fmoptions, 'Algorithm', 'sqp');
0149 Af = full(Af);
0150 Afeq = full(Afeq);
0151 otherwise
0152 error('fmincopf_solver: unknown algorithm specified in ''fmincon.alg'' option');
0153 end
0154 else
0155 fmoptions = optimset(fmoptions, 'LargeScale', 'off');
0156 Af = full(Af);
0157 Afeq = full(Afeq);
0158 end
0159
0160
0161
0162
0163 f_fcn = @(x)opf_costfcn(x, om);
0164 gh_fcn = @(x)opf_consfcn(x, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
0165 [x, f, info, Output, Lambda] = ...
0166 fmincon(f_fcn, x0, Af, bf, Afeq, bfeq, LB, UB, gh_fcn, fmoptions);
0167 success = (info > 0);
0168
0169
0170 Va = x(vv.i1.Va:vv.iN.Va);
0171 Vm = x(vv.i1.Vm:vv.iN.Vm);
0172 Pg = x(vv.i1.Pg:vv.iN.Pg);
0173 Qg = x(vv.i1.Qg:vv.iN.Qg);
0174 V = Vm .* exp(1j*Va);
0175
0176
0177
0178 bus(:, VA) = Va * 180/pi;
0179 bus(:, VM) = Vm;
0180 gen(:, PG) = Pg * baseMVA;
0181 gen(:, QG) = Qg * baseMVA;
0182 gen(:, VG) = Vm(gen(:, GEN_BUS));
0183
0184
0185 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);
0186 St = V(branch(:, T_BUS)) .* conj(Yt * V);
0187 branch(:, PF) = real(Sf) * baseMVA;
0188 branch(:, QF) = imag(Sf) * baseMVA;
0189 branch(:, PT) = real(St) * baseMVA;
0190 branch(:, QT) = imag(St) * baseMVA;
0191
0192
0193
0194 muSf = zeros(nl, 1);
0195 muSt = zeros(nl, 1);
0196 if ~isempty(il)
0197 muSf(il) = 2 * Lambda.ineqnonlin(1:nl2) .* branch(il, RATE_A) / baseMVA;
0198 muSt(il) = 2 * Lambda.ineqnonlin((1:nl2)+nl2) .* branch(il, RATE_A) / baseMVA;
0199 end
0200
0201
0202 bus(:, MU_VMAX) = Lambda.upper(vv.i1.Vm:vv.iN.Vm);
0203 bus(:, MU_VMIN) = Lambda.lower(vv.i1.Vm:vv.iN.Vm);
0204 gen(:, MU_PMAX) = Lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0205 gen(:, MU_PMIN) = Lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0206 gen(:, MU_QMAX) = Lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0207 gen(:, MU_QMIN) = Lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0208 bus(:, LAM_P) = Lambda.eqnonlin(nn.i1.Pmis:nn.iN.Pmis) / baseMVA;
0209 bus(:, LAM_Q) = Lambda.eqnonlin(nn.i1.Qmis:nn.iN.Qmis) / baseMVA;
0210 branch(:, MU_SF) = muSf / baseMVA;
0211 branch(:, MU_ST) = muSt / baseMVA;
0212
0213
0214 nlnN = getN(om, 'nln');
0215 nlt = length(ilt);
0216 ngt = length(igt);
0217 nbx = length(ibx);
0218
0219
0220 kl = find(Lambda.eqnonlin < 0);
0221 ku = find(Lambda.eqnonlin > 0);
0222 nl_mu_l = zeros(nlnN, 1);
0223 nl_mu_u = [zeros(2*nb, 1); muSf; muSt];
0224 nl_mu_l(kl) = -Lambda.eqnonlin(kl);
0225 nl_mu_u(ku) = Lambda.eqnonlin(ku);
0226
0227
0228 kl = find(Lambda.eqlin < 0);
0229 ku = find(Lambda.eqlin > 0);
0230
0231 mu_l = zeros(size(u));
0232 mu_l(ieq(kl)) = -Lambda.eqlin(kl);
0233 mu_l(igt) = Lambda.ineqlin(nlt+(1:ngt));
0234 mu_l(ibx) = Lambda.ineqlin(nlt+ngt+nbx+(1:nbx));
0235
0236 mu_u = zeros(size(u));
0237 mu_u(ieq(ku)) = Lambda.eqlin(ku);
0238 mu_u(ilt) = Lambda.ineqlin(1:nlt);
0239 mu_u(ibx) = Lambda.ineqlin(nlt+ngt+(1:nbx));
0240
0241 mu = struct( ...
0242 'var', struct('l', Lambda.lower, 'u', Lambda.upper), ...
0243 'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0244 'lin', struct('l', mu_l, 'u', mu_u) );
0245
0246 results = mpc;
0247 [results.bus, results.branch, results.gen, ...
0248 results.om, results.x, results.mu, results.f] = ...
0249 deal(bus, branch, gen, om, x, mu, f);
0250
0251 pimul = [ ...
0252 results.mu.nln.l - results.mu.nln.u;
0253 results.mu.lin.l - results.mu.lin.u;
0254 -ones(ny>0, 1);
0255 results.mu.var.l - results.mu.var.u;
0256 ];
0257 raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', Output);