Description of copf

0001 function [results, success, raw] = copf_solver(om, mpopt)
0002 %------------------------------  deprecated  ------------------------------
0003 %   OPF solver based on CONSTR to be removed in a future version.
0004 %--------------------------------------------------------------------------
0005 %COPF_SOLVER  Solves AC optimal power flow using CONSTR (Opt Tbx 1.x & 2.x).
0006 %
0007 %   [RESULTS, SUCCESS, RAW] = COPF_SOLVER(OM, MPOPT)
0008 %
0009 %   Inputs are an OPF model object and a MATPOWER options vector.
0010 %
0011 %   Outputs are a RESULTS struct, SUCCESS flag and RAW output struct.
0012 %
0013 %   RESULTS is a MATPOWER case struct (mpc) with the usual baseMVA, bus
0014 %   branch, gen, gencost fields, along with the following additional
0015 %   fields:
0016 %       .order      see 'help ext2int' for details of this field
0017 %       .x          final value of optimization variables (internal order)
0018 %       .f          final objective function value
0019 %       .mu         shadow prices on ...
0020 %           .var
0021 %               .l  lower bounds on variables
0022 %               .u  upper bounds on variables
0023 %           .nln
0024 %               .l  lower bounds on nonlinear constraints
0025 %               .u  upper bounds on nonlinear constraints
0026 %           .lin
0027 %               .l  lower bounds on linear constraints
0028 %               .u  upper bounds on linear constraints
0029 %
0030 %   SUCCESS     1 if solver converged successfully, 0 otherwise
0031 %
0032 %   RAW         raw output in form returned by MINOS
0033 %       .xr     final value of optimization variables
0034 %       .pimul  constraint multipliers
0035 %       .info   solver specific termination code
0036 %
0037 %   See also OPF, CONSTR, FUN_COPF, GRAD_COPF.
0038 
0039 %   MATPOWER
0040 %   $Id: copf_solver.m,v 1.22 2010/06/09 14:56:58 ray Exp $
0041 %   by Ray Zimmerman, PSERC Cornell
0042 %   and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0043 %   Copyright (c) 2000-2010 by Power System Engineering Research Center (PSERC)
0044 %
0045 %   This file is part of MATPOWER.
0046 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0047 %
0048 %   MATPOWER is free software: you can redistribute it and/or modify
0049 %   it under the terms of the GNU General Public License as published
0050 %   by the Free Software Foundation, either version 3 of the License,
0051 %   or (at your option) any later version.
0052 %
0053 %   MATPOWER is distributed in the hope that it will be useful,
0054 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0055 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0056 %   GNU General Public License for more details.
0057 %
0058 %   You should have received a copy of the GNU General Public License
0059 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0060 %
0061 %   Additional permission under GNU GPL version 3 section 7
0062 %
0063 %   If you modify MATPOWER, or any covered work, to interface with
0064 %   other modules (such as MATLAB code and MEX-files) available in a
0065 %   MATLAB(R) or comparable environment containing parts covered
0066 %   under other licensing terms, the licensors of MATPOWER grant
0067 %   you additional permission to convey the resulting work.
0068 
0069 %%----- initialization -----
0070 %% define named indices into data matrices
0071 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0072     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0073 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0074     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0075     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0076 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0077     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0078     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0079 
0080 %% options
0081 verbose = mpopt(31);    %% VERBOSE
0082 
0083 %% unpack data
0084 mpc = get_mpc(om);
0085 [baseMVA, bus, gen, branch] = ...
0086     deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch);
0087 vv = get_idx(om);
0088 
0089 %% problem dimensions
0090 nb = size(bus, 1);          %% number of buses
0091 nl = size(branch, 1);       %% number of branches
0092 ny = getN(om, 'var', 'y');  %% number of piece-wise linear costs
0093 
0094 %% bounds on optimization vars
0095 [x0, LB, UB] = getv(om);
0096 
0097 %% linear constraints
0098 %% WORKAROUND: Add bounds on all vars to A, l, u
0099 nxyz = length(x0);
0100 om2 = om;
0101 om2 = add_constraints(om2, 'varlims', speye(nxyz, nxyz), LB, UB);
0102 [vv, ll, nn] = get_idx(om2);
0103 [A, l, u] = linear_constraints(om2);
0104 
0105 %% so, can we do anything good about lambda initialization?
0106 if all(bus(:, LAM_P) == 0)
0107   bus(:, LAM_P) = (10)*ones(nb, 1);
0108 end
0109 
0110 %% split l <= A*x <= u into less than, equal to, greater than, and
0111 %% doubly-bounded sets
0112 ieq = find( abs(u-l) <= eps );          %% equality
0113 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0114 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0115 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0116 Af  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0117 bf  = [ u(ilt);   -l(igt);     u(ibx);    -l(ibx)];
0118 Afeq = A(ieq, :);
0119 bfeq = u(ieq);
0120 
0121 %% build admittance matrices
0122 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0123 
0124 %% find branches with flow limits
0125 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0126 nl2 = length(il);           %% number of constrained lines
0127 
0128 %% set mpopt defaults
0129 mpopt(15)  = 2 * nb + length(bfeq); %% set number of equality constraints
0130 if mpopt(19) == 0           %% CONSTR_MAX_IT
0131   mpopt(19) = 150 + 2*nb;
0132 end
0133 
0134 %% set up options for Optim Tbx's constr
0135 otopt = foptions;           %% get default options for constr
0136 otopt(1) = (verbose > 0);   %% set verbose flag appropriately
0137 % otopt(9) = 1;             %% check user supplied gradients?
0138 otopt(2)  = mpopt(17);      %% termination tolerance on 'x'
0139 otopt(3)  = mpopt(18);      %% termination tolerance on 'F'
0140 otopt(4)  = mpopt(16);      %% termination tolerance on constraint violation
0141 otopt(13) = mpopt(15);      %% number of equality constraints
0142 otopt(14) = mpopt(19);      %% maximum number of iterations
0143 
0144 %% run optimization
0145 [x, otopt, lambda] = constr('fun_copf', x0, otopt, [], [], 'grad_copf', ...
0146         om2, Ybus, Yf(il,:), Yt(il,:), Afeq, bfeq, Af, bf, mpopt, il);
0147 %% get final objective function value & constraint values
0148 [f, g] = feval('fun_copf', x, om2, Ybus, Yf(il,:), Yt(il,:), Afeq, bfeq, Af, bf, mpopt, il);
0149 
0150 %% check for convergence
0151 if otopt(10) >= otopt(14) || max(abs(g(1:otopt(13)))) > otopt(4) ...
0152                    || max(g((otopt(13)+1):length(g))) > otopt(4)
0153   success = 0;        %% did NOT converge
0154 else
0155   success = 1;        %% DID converge
0156 end
0157 info = success;
0158 
0159 %% update solution data
0160 Va = x(vv.i1.Va:vv.iN.Va);
0161 Vm = x(vv.i1.Vm:vv.iN.Vm);
0162 Pg = x(vv.i1.Pg:vv.iN.Pg);
0163 Qg = x(vv.i1.Qg:vv.iN.Qg);
0164 V = Vm .* exp(1j*Va);
0165 
0166 %%-----  calculate return values  -----
0167 %% update voltages & generator outputs
0168 bus(:, VA) = Va * 180/pi;
0169 bus(:, VM) = Vm;
0170 gen(:, PG) = Pg * baseMVA;
0171 gen(:, QG) = Qg * baseMVA;
0172 gen(:, VG) = Vm(gen(:, GEN_BUS));
0173 
0174 %% compute branch flows
0175 Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
0176 St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
0177 branch(:, PF) = real(Sf) * baseMVA;
0178 branch(:, QF) = imag(Sf) * baseMVA;
0179 branch(:, PT) = real(St) * baseMVA;
0180 branch(:, QT) = imag(St) * baseMVA;
0181 
0182 %% package up results
0183 nA = length(u);
0184 neq = length(ieq);
0185 nlt = length(ilt);
0186 ngt = length(igt);
0187 nbx = length(ibx);
0188 
0189 %% extract multipliers (lambda is ordered as Pmis, Qmis, Afeq, Sf, St, Af)
0190 %% nonlinear constraints
0191 inln = [(1:2*nb), (1:2*nl2) + 2*nb+neq];
0192 kl = find(lambda(inln) < 0);
0193 ku = find(lambda(inln) > 0);
0194 nl_mu_l = zeros(2*(nb+nl2), 1);
0195 nl_mu_u = zeros(2*(nb+nl2), 1);
0196 nl_mu_l(kl) = -lambda(inln(kl));
0197 nl_mu_u(ku) =  lambda(inln(ku));
0198 
0199 %% linear constraints
0200 ilin = [(1:neq)+2*nb, (1:(nlt+ngt+2*nbx)) + 2*nb+neq+2*nl2];
0201 kl = find(lambda(ilin(1:neq)) < 0);
0202 ku = find(lambda(ilin(1:neq)) > 0);
0203 
0204 mu_l = zeros(nA, 1);
0205 mu_l(ieq) = -lambda(ilin(1:neq));
0206 mu_l(ieq(ku)) = 0;
0207 mu_l(igt) = lambda(ilin(neq+nlt+(1:ngt)));
0208 mu_l(ibx) = lambda(ilin(neq+nlt+ngt+nbx+(1:nbx)));
0209 
0210 mu_u = zeros(nA, 1);
0211 mu_u(ieq) = lambda(ilin(1:neq));
0212 mu_u(ieq(kl)) = 0;
0213 mu_u(ilt) = lambda(ilin(neq+(1:nlt)));
0214 mu_u(ibx) = lambda(ilin(neq+nlt+ngt+(1:nbx)));
0215 
0216 %% variable bounds
0217 muLB = mu_l(ll.i1.varlims:ll.iN.varlims);
0218 muUB = mu_u(ll.i1.varlims:ll.iN.varlims);
0219 mu_l(ll.i1.varlims:ll.iN.varlims) = [];
0220 mu_u(ll.i1.varlims:ll.iN.varlims) = [];
0221 
0222 %% line constraint is actually on square of limit
0223 %% so we must fix multipliers
0224 muSf = zeros(nl, 1);
0225 muSt = zeros(nl, 1);
0226 muSf(il) = 2 * nl_mu_u((1:nl2)+2*nb    ) .* branch(il, RATE_A) / baseMVA;
0227 muSt(il) = 2 * nl_mu_u((1:nl2)+2*nb+nl2) .* branch(il, RATE_A) / baseMVA;
0228 
0229 %% resize mu for nonlinear constraints
0230 nl_mu_l = [nl_mu_l(1:2*nb); zeros(2*nl, 1)];
0231 nl_mu_u = [nl_mu_u(1:2*nb); muSf; muSt];
0232 
0233 %% update Lagrange multipliers
0234 bus(:, MU_VMAX)  = muUB(vv.i1.Vm:vv.iN.Vm);
0235 bus(:, MU_VMIN)  = muLB(vv.i1.Vm:vv.iN.Vm);
0236 gen(:, MU_PMAX)  = muUB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0237 gen(:, MU_PMIN)  = muLB(vv.i1.Pg:vv.iN.Pg) / baseMVA;
0238 gen(:, MU_QMAX)  = muUB(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0239 gen(:, MU_QMIN)  = muLB(vv.i1.Qg:vv.iN.Qg) / baseMVA;
0240 bus(:, LAM_P)    = (nl_mu_u(nn.i1.Pmis:nn.iN.Pmis) - nl_mu_l(nn.i1.Pmis:nn.iN.Pmis)) / baseMVA;
0241 bus(:, LAM_Q)    = (nl_mu_u(nn.i1.Qmis:nn.iN.Qmis) - nl_mu_l(nn.i1.Qmis:nn.iN.Qmis)) / baseMVA;
0242 branch(:, MU_SF) = muSf / baseMVA;
0243 branch(:, MU_ST) = muSt / baseMVA;
0244 
0245 mu = struct( ...
0246   'var', struct('l', muLB, 'u', muUB), ...
0247   'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ...
0248   'lin', struct('l', mu_l, 'u', mu_u) );
0249 
0250 results = mpc;
0251 [results.bus, results.branch, results.gen, ...
0252     results.om, results.x, results.mu, results.f] = ...
0253         deal(bus, branch, gen, om, x, mu, f);
0254 
0255 pimul = [ ...
0256   results.mu.nln.l - results.mu.nln.u;
0257   results.mu.lin.l - results.mu.lin.u;
0258   -ones(ny>0, 1);
0259   results.mu.var.l - results.mu.var.u;
0260 ];
0261 raw = struct('xr', x, 'pimul', pimul, 'info', info);
copf_solver

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
