consfmin.m

该程序是计算最优潮流的matlab工具箱。可以很好的求解目标函数不同的最优潮流问题。

function [g, geq, dg, dgeq] = consfmin(x, baseMVA, bus, gen, gencost, branch, areas, Ybus, Yf, Yt, mpopt, parms, ccost)
%CONSFMIN  Evaluates nonlinear constraints and their Jacobian for OPF.
%   [g, geq, dg, dgeq] = consfmin(x, baseMVA, bus, gen, gencost, ...
%                                  branch, areas, Ybus, Yf, Yt, mpopt)

%   MATPOWER
%   $Id: consfmin.m,v 1.4 2004/09/07 18:27:37 ray Exp $
%   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
%   and Ray Zimmerman, PSERC Cornell
%   Copyright (c) 1996-2004 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/matpower/ for more info.

%%----- initialize -----
%% define named indices into data matrices
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ...
  GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = 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] = idx_brch;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, N, COST] = idx_cost;

%% constant
j = sqrt(-1);

%% generator info
on = find(gen(:, GEN_STATUS) > 0);    %% which generators are on?
gbus = gen(on, GEN_BUS);        %% what buses are they at?

%% sizes of things
nb = size(bus, 1);
nl = size(branch, 1);
%npv  = length(pv);
%npq  = length(pq);
ng = length(on);            %% number of generators that are turned on
br = find(branch(:, BR_STATUS));
nbr = length(br);

%% set up indexing for x
j1 = 1;      j4 = nb;             %% j1:j4 - V angle of all buses
j5 = j4 + 1;  j6  = j4 + nb;      %% j5:j6  - V mag of all buses
j7 = j6 + 1;  j8  = j6 + ng;      %% j7:j8  - P of generators
j9 = j8 + 1;  j10  = j8 + ng;      %% j9:j10  - Q of generators

nb = parms(1);
ng = parms(2);
nl = parms(3);
ny = parms(4);
nx = parms(5);
nvl = parms(6);
nz = parms(7);
nxyz = parms(8);
thbas = parms(9);
thend = parms(10);
vbas = parms(11);
vend = parms(12);
pgbas = parms(13);
pgend = parms(14);
qgbas = parms(15);
qgend = parms(16);
ybas = parms(17);
yend = parms(18);
zbas = parms(19);
zend = parms(20);
pmsmbas = parms(21);
pmsmend = parms(22);
qmsmbas = parms(23);
qmsmend = parms(24);
sfbas = parms(25);
sfend = parms(26);
stbas = parms(27);
stend = parms(28);

%% grab Pg & Qg
Pg = x(j7:j8);                %% active generation in p.u.
Qg = x(j9:j10);                %% reactive generation in p.u.

%% put Pg & Qg back in gen
gen(on, PG) = Pg * baseMVA;          %% active generation in MW
gen(on, QG) = Qg * baseMVA;          %% reactive generation in MVAr
 
%% rebuild Sbus
Sbus = makeSbus(baseMVA, bus, gen);    %% net injected power in p.u.

%% ----- evaluate constraints -----
%% reconstruct V
Va = zeros(nb, 1);
Va = x(j1:j4);
Vm = x(j5:j6);
V = Vm .* exp(j * Va);

%% evaluate power flow equations
mis = V .* conj(Ybus * V) - Sbus;

%% compute branch power flows
Sf = V(branch(br, F_BUS)) .* conj(Yf(br, :) * V);  %% complex power injected at "from" bus (p.u.)
St = V(branch(br, T_BUS)) .* conj(Yt(br, :) * V);  %% complex power injected at "to" bus (p.u.)

%% compute constraint function values

%% First the equality constraints (power flow)
geq = [ real(mis);              %% active power mismatch for all buses
        imag(mis) ];            %% reactive power mismatch for all buses

%% Then the inequality constraints (branch limits)
if mpopt(24) == 1   %% P limit (Pan Wei)
  g = [ real(Sf) - branch(br, RATE_A)/baseMVA;  %% branch apparent power limits (from bus)
        real(St) - branch(br, RATE_A)/baseMVA ]; %% branch apparent power limits (to bus)
else                %% |S| limit
  g = [ abs(Sf) - branch(br, RATE_A)/baseMVA;  %% branch apparent power limits (from bus)
        abs(St) - branch(br, RATE_A)/baseMVA ]; %% branch apparent power limits (to bus)
end

%%----- evaluate partials of constraints -----
if nargout > 2
  %% compute partials of injected bus powers
  [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V);    %% w.r.t. V
  dSbus_dPg = sparse(gbus, 1:ng, -1, nb, ng);    %% w.r.t. Pg
  dSbus_dQg = sparse(gbus, 1:ng, -j, nb, ng);    %% w.r.t. Qg
  
  %% compute partials of line flows w.r.t. V
  [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V);

  %% line limits are w.r.t apparent power, so compute partials of apparent power
  [dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm] = ...
            dAbr_dV(dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St);

  %% construct Jacobian of equality constraints (load flows) and transpose it
  dgeq = [
    %% equality constraints
    real(dSbus_dVa), real(dSbus_dVm), ...
          real(dSbus_dPg), real(dSbus_dQg), sparse(nb, ny+nz);   %% P mismatch
    imag(dSbus_dVa), imag(dSbus_dVm), ...
          imag(dSbus_dPg), imag(dSbus_dQg), sparse(nb, ny+nz);   %% Q mismatch
   ]'; 
  %% construct Jacobian of inequality constraints (branch limits)
  %% and transpose it so fmincon likes it
  if mpopt(24) == 1     %% P limit (Pan Wei)
    dg = [
      real(dSf_dVa(br,:)), real(dSf_dVm(br,:)), sparse(nbr,2*ng+ny+nz);		%% Pf limit
      real(dSt_dVa(br,:)), real(dSt_dVm(br,:)), sparse(nbr,2*ng+ny+nz);		%% Pt limit
    ]';
  else                  %% |S| limit
    dg = [
      dAf_dVa(br,:), dAf_dVm(br,:), sparse(nbr,2*ng+ny+nz);  %% |Sf| limit
      dAt_dVa(br,:), dAt_dVm(br,:), sparse(nbr,2*ng+ny+nz);  %% |St| limit
    ]';
  end  

  %% the following lines can be removed if/when fmincon
  %% supports sparse matrices for non-linearly constrained problems
  dgeq = full(dgeq);
  dg = full(dg);
end

return;