consfmin.m

用于潮流和最优潮流计算的最新版本 对于初学者很有用

function [g, geq, dg, dgeq] = consfmin(x, baseMVA, bus, gen, gencost, branch, areas, Ybus, Yf, Yt, mpopt, parms, ccost, N, fparm, H, Cw)
%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, parms, ccost, N, fparm, H, Cw)

%   MATPOWER
%   $Id: consfmin.m,v 1.12 2007/06/26 16:05:21 ray Exp $
%   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
%   and Ray Zimmerman, PSERC Cornell
%   Copyright (c) 1996-2006 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, 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;

%% constant
j = sqrt(-1);

%% unpack needed parameters
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(pgbas:pgend);            %% active generation in p.u.
Qg = x(qgbas:qgend);            %% reactive generation in p.u.

%% put Pg & Qg back in gen
gen(:, PG) = Pg * baseMVA;      %% active generation in MW
gen(:, 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(thbas:thend);
Vm = x(vbas:vend);
V = Vm .* exp(j * Va);

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

%%----- evaluate 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 flow limits)
if mpopt(24) == 2       %% current magnitude limit, |I|
  g = [ abs(Yf*V) - branch(:, RATE_A)/baseMVA;    %% branch current limits (from bus)
        abs(Yt*V) - branch(:, RATE_A)/baseMVA ];  %% branch current limits (to bus)
else
    %% compute branch power flows
    Sf = V(branch(:, F_BUS)) .* conj(Yf * V);   %% complex power injected at "from" bus (p.u.)
    St = V(branch(:, T_BUS)) .* conj(Yt * V);   %% complex power injected at "to" bus (p.u.)
    if mpopt(24) == 1   %% active power limit, P (Pan Wei)
      g = [ real(Sf) - branch(:, RATE_A)/baseMVA;   %% branch real power limits (from bus)
            real(St) - branch(:, RATE_A)/baseMVA ]; %% branch real power limits (to bus)
    else                %% apparent power limit, |S|
      g = [ abs(Sf) - branch(:, RATE_A)/baseMVA;    %% branch apparent power limits (from bus)
            abs(St) - branch(:, RATE_A)/baseMVA ];  %% branch apparent power limits (to bus)
    end
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(gen(:, GEN_BUS), 1:ng, -1, nb, ng);    %% w.r.t. Pg
  dSbus_dQg = sparse(gen(:, GEN_BUS), 1:ng, -j, nb, ng);    %% w.r.t. Qg
  
  %% construct Jacobian of equality constraints (power flow) 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
   ]'; 

  %% compute partials of Flows w.r.t. V
  if mpopt(24) == 2     %% current
    [dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft] = dIbr_dV(branch, Yf, Yt, V);
  else                  %% power
    [dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft] = dSbr_dV(branch, Yf, Yt, V);
  end
  if mpopt(24) == 1     %% real part of flow (active power)
    df_dVa = real(dFf_dVa);
    df_dVm = real(dFf_dVm);
    dt_dVa = real(dFt_dVa);
    dt_dVm = real(dFt_dVm);
  else                  %% magnitude of flow (of complex power or current)
    [df_dVa, df_dVm, dt_dVa, dt_dVm] = ...
            dAbr_dV(dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft);
  end

  %% construct Jacobian of inequality constraints (branch limits)
  %% and transpose it so fmincon likes it
  dg = [
    df_dVa, df_dVm, sparse(nl,2*ng+ny+nz);  %% "from" flow limit
    dt_dVa, dt_dVm, sparse(nl,2*ng+ny+nz);  %% "to" flow limit
  ]';

  %% 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;