www.gusucode.com > 基于lingo求所以解,对潮流计算求出所有解 > matpower4.1/mips6opf_solver.m
function [results, success, raw] = mips6opf_solver(om, mpopt) %------------------------------ deprecated ------------------------------ % MATLAB 6.x support to be removed in a future version. %-------------------------------------------------------------------------- %MIPS6OPF_SOLVER Solves AC optimal power flow using MIPS (for MATLAB 6.x). % % [RESULTS, SUCCESS, RAW] = MIPS6OPF_SOLVER(OM, MPOPT) % % Inputs are an OPF model object and a MATPOWER options vector. % % Outputs are a RESULTS struct, SUCCESS flag and RAW output struct. % % RESULTS is a MATPOWER case struct (mpc) with the usual baseMVA, bus % branch, gen, gencost fields, along with the following additional % fields: % .order see 'help ext2int' for details of this field % .x final value of optimization variables (internal order) % .f final objective function value % .mu shadow prices on ... % .var % .l lower bounds on variables % .u upper bounds on variables % .nln % .l lower bounds on nonlinear constraints % .u upper bounds on nonlinear constraints % .lin % .l lower bounds on linear constraints % .u upper bounds on linear constraints % % SUCCESS 1 if solver converged successfully, 0 otherwise % % RAW raw output in form returned by MINOS % .xr final value of optimization variables % .pimul constraint multipliers % .info solver specific termination code % .output solver specific output information % % See also OPF, MIPS6. % MATPOWER % $Id: mips6opf_solver.m,v 1.16 2010/06/09 14:56:58 ray Exp $ % by Ray Zimmerman, PSERC Cornell % and Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales % Copyright (c) 2000-2010 by Power System Engineering Research Center (PSERC) % % This file is part of MATPOWER. % See http://www.pserc.cornell.edu/matpower/ for more info. % % MATPOWER is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published % by the Free Software Foundation, either version 3 of the License, % or (at your option) any later version. % % MATPOWER is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. % % Additional permission under GNU GPL version 3 section 7 % % If you modify MATPOWER, or any covered work, to interface with % other modules (such as MATLAB code and MEX-files) available in a % MATLAB(R) or comparable environment containing parts covered % under other licensing terms, the licensors of MATPOWER grant % you additional permission to convey the resulting work. %%----- initialization ----- %% define named indices into data matrices [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost; %% options verbose = mpopt(31); %% VERBOSE feastol = mpopt(81); %% PDIPM_FEASTOL gradtol = mpopt(82); %% PDIPM_GRADTOL comptol = mpopt(83); %% PDIPM_COMPTOL costtol = mpopt(84); %% PDIPM_COSTTOL max_it = mpopt(85); %% PDIPM_MAX_IT max_red = mpopt(86); %% SCPDIPM_RED_IT step_control = (mpopt(11) == 565); %% OPF_ALG == 565, MIPS-sc if feastol == 0 feastol = mpopt(16); %% = OPF_VIOLATION by default end opt = struct( 'feastol', feastol, ... 'gradtol', gradtol, ... 'comptol', comptol, ... 'costtol', costtol, ... 'max_it', max_it, ... 'max_red', max_red, ... 'step_control', step_control, ... 'cost_mult', 1e-4, ... 'verbose', verbose ); %% unpack data mpc = get_mpc(om); [baseMVA, bus, gen, branch, gencost] = ... deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch, mpc.gencost); [vv, ll, nn] = get_idx(om); %% problem dimensions nb = size(bus, 1); %% number of buses nl = size(branch, 1); %% number of branches ny = getN(om, 'var', 'y'); %% number of piece-wise linear costs %% linear constraints [A, l, u] = linear_constraints(om); %% bounds on optimization vars [x0, xmin, xmax] = getv(om); %% build admittance matrices [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch); %% try to select an interior initial point ll = xmin; uu = xmax; ll(xmin == -Inf) = -1e10; %% replace Inf with numerical proxies uu(xmax == Inf) = 1e10; x0 = (ll + uu) / 2; Varefs = bus(bus(:, BUS_TYPE) == REF, VA) * (pi/180); x0(vv.i1.Va:vv.iN.Va) = Varefs(1); %% angles set to first reference angle if ny > 0 ipwl = find(gencost(:, MODEL) == PW_LINEAR); % PQ = [gen(:, PMAX); gen(:, QMAX)]; % c = totcost(gencost(ipwl, :), PQ(ipwl)); c = gencost(sub2ind(size(gencost), ipwl, NCOST+2*gencost(ipwl, NCOST))); %% largest y-value in CCV data x0(vv.i1.y:vv.iN.y) = max(c) + 0.1 * abs(max(c)); % x0(vv.i1.y:vv.iN.y) = c + 0.1 * abs(c); end %% find branches with flow limits il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10); nl2 = length(il); %% number of constrained lines %%----- run opf ----- f_fcn = @opf_costfcn; gh_fcn = @opf_consfcn; ipm_hessian = @opf_hessfcn; [x, f, info, Output, Lambda] = ... mips6(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, ipm_hessian, opt, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il); success = (info > 0); %% update solution data Va = x(vv.i1.Va:vv.iN.Va); Vm = x(vv.i1.Vm:vv.iN.Vm); Pg = x(vv.i1.Pg:vv.iN.Pg); Qg = x(vv.i1.Qg:vv.iN.Qg); V = Vm .* exp(1j*Va); %%----- calculate return values ----- %% update voltages & generator outputs bus(:, VA) = Va * 180/pi; bus(:, VM) = Vm; gen(:, PG) = Pg * baseMVA; gen(:, QG) = Qg * baseMVA; gen(:, VG) = Vm(gen(:, GEN_BUS)); %% compute branch flows Sf = V(branch(:, F_BUS)) .* conj(Yf * V); %% cplx pwr at "from" bus, p.u. St = V(branch(:, T_BUS)) .* conj(Yt * V); %% cplx pwr at "to" bus, p.u. branch(:, PF) = real(Sf) * baseMVA; branch(:, QF) = imag(Sf) * baseMVA; branch(:, PT) = real(St) * baseMVA; branch(:, QT) = imag(St) * baseMVA; %% line constraint is actually on square of limit %% so we must fix multipliers muSf = zeros(nl, 1); muSt = zeros(nl, 1); if ~isempty(il) muSf(il) = 2 * Lambda.ineqnonlin(1:nl2) .* branch(il, RATE_A) / baseMVA; muSt(il) = 2 * Lambda.ineqnonlin((1:nl2)+nl2) .* branch(il, RATE_A) / baseMVA; end %% update Lagrange multipliers bus(:, MU_VMAX) = Lambda.upper(vv.i1.Vm:vv.iN.Vm); bus(:, MU_VMIN) = Lambda.lower(vv.i1.Vm:vv.iN.Vm); gen(:, MU_PMAX) = Lambda.upper(vv.i1.Pg:vv.iN.Pg) / baseMVA; gen(:, MU_PMIN) = Lambda.lower(vv.i1.Pg:vv.iN.Pg) / baseMVA; gen(:, MU_QMAX) = Lambda.upper(vv.i1.Qg:vv.iN.Qg) / baseMVA; gen(:, MU_QMIN) = Lambda.lower(vv.i1.Qg:vv.iN.Qg) / baseMVA; bus(:, LAM_P) = Lambda.eqnonlin(nn.i1.Pmis:nn.iN.Pmis) / baseMVA; bus(:, LAM_Q) = Lambda.eqnonlin(nn.i1.Qmis:nn.iN.Qmis) / baseMVA; branch(:, MU_SF) = muSf / baseMVA; branch(:, MU_ST) = muSt / baseMVA; %% package up results nlnN = getN(om, 'nln'); %% extract multipliers for nonlinear constraints kl = find(Lambda.eqnonlin < 0); ku = find(Lambda.eqnonlin > 0); nl_mu_l = zeros(nlnN, 1); nl_mu_u = [zeros(2*nb, 1); muSf; muSt]; nl_mu_l(kl) = -Lambda.eqnonlin(kl); nl_mu_u(ku) = Lambda.eqnonlin(ku); mu = struct( ... 'var', struct('l', Lambda.lower, 'u', Lambda.upper), ... 'nln', struct('l', nl_mu_l, 'u', nl_mu_u), ... 'lin', struct('l', Lambda.mu_l, 'u', Lambda.mu_u) ); results = mpc; [results.bus, results.branch, results.gen, ... results.om, results.x, results.mu, results.f] = ... deal(bus, branch, gen, om, x, mu, f); pimul = [ ... results.mu.nln.l - results.mu.nln.u; results.mu.lin.l - results.mu.lin.u; -ones(ny>0, 1); results.mu.var.l - results.mu.var.u; ]; raw = struct('xr', x, 'pimul', pimul, 'info', info, 'output', Output);