equ and ine constraints.m

求最优潮流的matlab程序

[Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);

%% 根据传进来的x值进行约束计算
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)
    
    
    
    
%%%——————-直流潮流的约束形式表达————————
% gen 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);
ng = length(on);                        %% number of generators that are turned on

%% initial state
Va  = bus(:, VA) * (pi/180);
Pg  = gen(on, PG) / baseMVA;

%% get bus index lists of each type of bus
[ref, pv, pq] = bustypes(bus, gen);

%% build B matrices and phase shift injections
[B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch);  %%在这里把计算的潮流等式和支路潮流计算的不等式也写出来

AA = [
        sparse(1, ref, 1, 1, nb+ng+nc);                     %% reference angle
        B,  -sparse(gen(on, GEN_BUS), 1:ng, ones(ng, 1), nb, ng), sparse(nb, nc);
                                                            %% real power flow eqns
        sparse(ng, nb), -speye(ng, ng), sparse(ng, nc);     %% lower limit on Pg
        sparse(ng, nb), speye(ng, ng), sparse(ng, nc);      %% upper limit on Pg
        Bf, sparse(nl, ng+nc);                              %% flow limit on Pf
        -Bf, sparse(nl, ng+nc);                             %% flow limit on Pt
        Acc;                                                %% cost constraints
];

bb = [
        Va(ref);                                            %% reference angle
        -(bus(:, PD) + bus(:, GS)) / baseMVA - Pbusinj;     %% real power flow eqns
        -gen(on, PMIN) / baseMVA;                           %% lower limit on Pg
        gen(on, PMAX) / baseMVA;                            %% upper limit on Pg
        branch(:, RATE_A) / baseMVA - Pfinj;                %% flow limit on Pf
        branch(:, RATE_A) / baseMVA + Pfinj;                %% flow limit on Pt
        bcc;                                                %% cost constraints
];


%%%——————————线性规划约束表达————
    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 line flow constraints
    br = find(branch(:, BR_STATUS));

        flow_limit = [
            abs(Yf(br, :)*V) - branch(br, RATE_A) / baseMVA;    %% from bus
            abs(Yt(br, :)*V) - branch(br, RATE_A) / baseMVA;    %% to bus
         ];
         
    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;
    ];