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