fun_std.m

matlab语言 编写的电力系统仿真代码

function [f, g] = fun_std(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt)
%FUN_STD  Evaluates objective function & constraints for OPF.
%   [f, g] = fun_std(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt)

%   MATPOWER
%   $Id: fun_std.m,v 1.6 2004/09/08 12:37:36 ray Exp $
%   by 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 gen, branch 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] = 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?

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

%% set up indexing for x
j1 = 1;         j2  = npv;              %% j1:j2    - V angle of pv buses
j3 = j2 + 1;    j4  = j2 + npq;         %% j3:j4    - V angle of pq 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

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

%%----- evaluate objective function -----
%% put Pg & Qg back in gen
gen(on, PG) = Pg * baseMVA;                 %% active generation in MW
gen(on, QG) = Qg * baseMVA;                 %% reactive generation in MVAr

%% compute objective value
[pcost, qcost] = pqcost(gencost, size(gen, 1), on);
f = sum(    [totcost(pcost, gen(on, PG)); ...   %% cost of Pg
             totcost(qcost, gen(on, QG)) ]  );  %% cost of Qg, empty if no qcost

%%----- evaluate constraints -----
if nargout > 1
    %% reconstruct V
    Va = zeros(nb, 1);
    Va([ref; pv; pq]) = [angle(V(ref)); x(j1:j2); x(j3:j4)];
    Vm = x(j5:j6);
    V = Vm .* exp(j * Va);
    
    %% rebuild Sbus
    Sbus = makeSbus(baseMVA, bus, gen);     %% net injected power in p.u.
    
    %% evaluate power flow equations
    mis = V .* conj(Ybus * V) - Sbus;
    
    %% compute branch power flows
    br = find(branch(:, BR_STATUS));
    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 line flow constraints
    if mpopt(24) == 1   %% branch active power limits
        flow_limit = [
            real(Sf) - branch(br, RATE_A) / baseMVA;    %% from bus
            real(St) - branch(br, RATE_A) / baseMVA;    %% to bus
        ];
    else                %% branch apparent power limits
        flow_limit = [
            abs(Sf) - branch(br, RATE_A) / baseMVA;     %% from bus
            abs(St) - branch(br, RATE_A) / baseMVA;     %% to bus
        ];
    end
    
    %% compute constraint function values
    g = [
        %% equality constraints
        real(mis);                          %% active power mismatch for all buses
        imag(mis);                          %% reactive power mismatch for all buses

        %% inequality constraints (variable limits, voltage & generation)
        bus(:, VMIN) - Vm;                  %% lower voltage limit for var V
        Vm - bus(:, VMAX);                  %% upper voltage limit for var V
        gen(on, PMIN) / baseMVA - Pg;       %% lower generator P limit
        Pg - gen(on, PMAX) / baseMVA;       %% upper generator P limit
        gen(on, QMIN) / baseMVA - Qg;       %% lower generator Q limit
        Qg - gen(on, QMAX) / baseMVA;       %% upper generator Q limit

        %% inequality constraints (line flow limits)
        flow_limit;
    ];
end

return;