fm_mass.m

电力系统的psat

function fm_mass(flag)
% FM_MASS define dynamic model for synchronous generator shaft
%
% FM_MASS(FLAG)
%       FLAG -> 0 initialization
%       FLAG -> 3 differential equations
%       FLAG -> 4 state matrix
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Version:   1.0.0
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2005 Federico Milano
%
% This toolbox 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 2.0 of the License, or
% (at your option) any later version.
%
% This toolbox is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANDABILITY 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 this toolbox; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307,
% USA.

global Mass DAE Syn Settings

delta = DAE.x(Syn.delta(Mass.syn));
omega = DAE.x(Syn.omega(Mass.syn));
dHP = DAE.x(Mass.delta_HP);
oHP = DAE.x(Mass.omega_HP);
dIP = DAE.x(Mass.delta_IP);
oIP = DAE.x(Mass.omega_IP);
dLP = DAE.x(Mass.delta_LP);
oLP = DAE.x(Mass.omega_LP);
dEX = DAE.x(Mass.delta_EX);
oEX = DAE.x(Mass.omega_EX);

Tm = Syn.pm(Mass.syn);
iM = 1./Syn.con(Mass.syn,18);

iMhp = 1./Mass.con(:,2);
iMip = 1./Mass.con(:,3);
iMlp = 1./Mass.con(:,4);
iMex = 1./Mass.con(:,5);

Dhp = Mass.con(:,6);
Dip = Mass.con(:,7);
Dlp = Mass.con(:,8);
Dex = Mass.con(:,9);

D12 = Mass.con(:,10);
D23 = Mass.con(:,11);
D34 = Mass.con(:,12);
D45 = Mass.con(:,13);

K12 = Mass.con(:,14);
K23 = Mass.con(:,15);
K34 = Mass.con(:,16);
K45 = Mass.con(:,17);

Wn = Settings.rad;

switch flag
case 0

    DAE.x(Mass.omega_HP) = omega;
    DAE.x(Mass.omega_IP) = omega;
    DAE.x(Mass.omega_LP) = omega;
    DAE.x(Mass.omega_EX) = omega;

    for i = 1:Mass.n
        DAE.x(Mass.delta_EX(i)) = delta(i);
        DAE.x(Mass.delta_LP(i)) = (Tm(i)+K34(i)*delta(i))/K34(i);
        DAE.x(Mass.delta_IP(i)) = (Tm(i)+K23(i)*DAE.x(Mass.delta_LP(i)))/K23(i);
        DAE.x(Mass.delta_HP(i)) = (Tm(i)+K12(i)*DAE.x(Mass.delta_IP(i)))/K12(i);
        fm_disp(['Initialization of Shaft Speeds & Angles  ( -> Bus',num2str(Syn.bus(Mass.syn(i))),') completed.'],1)
    end

case 3

    DAE.f(Mass.delta_HP) = Wn*(oHP-1);
    DAE.f(Mass.omega_HP) = (Tm-Dhp.*(oHP-1)-D12.*(oHP-oIP)-K12.*dHP+K12.*dIP).*iMhp;
    DAE.f(Mass.delta_IP) = Wn*(oIP-1);
    DAE.f(Mass.omega_IP) = (-Dip.*(oIP-1)-D12.*(oIP-oHP)-D23.*(oIP-oLP)+K12.*dHP-(K12+K23).*dIP+K23.*dLP).*iMip;
    DAE.f(Mass.delta_LP) = Wn*(oLP-1);
    DAE.f(Mass.omega_LP) = (-Dlp.*(oLP-1)-D23.*(oLP-oIP)-D34.*(oLP-omega)+K23.*dIP-(K23+K34).*dLP+K34.*delta).*iMlp;
    DAE.f(Syn.omega(Mass.syn)) = DAE.f(Syn.omega(Mass.syn)) + (-D34.*(omega-oLP)-D45.*(omega-oEX)+K34.*dLP -(K34+K45).*delta+K45.*dEX-Tm).*iM;
    DAE.f(Mass.delta_EX) = Wn*(oEX-1);
    DAE.f(Mass.omega_EX) = (-Dex.*(oEX-1)-D45.*(oEX-omega)+K45.*delta-K45.*dEX).*iMex;

case 4

    for i = 1:Mass.n
        k = Mass.syn(i);
        DAE.Fx(Mass.delta_HP(i),Mass.omega_HP(i)) = Wn;
        DAE.Fx(Mass.omega_HP(i),Mass.delta_HP(i)) = -K12(i)*iMhp(i);
        DAE.Fx(Mass.omega_HP(i),Mass.omega_HP(i)) = -(Dhp(i)+D12(i))*iMhp(i);
        DAE.Fx(Mass.omega_HP(i),Mass.delta_IP(i)) = K12(i)*iMhp(i);
        DAE.Fx(Mass.delta_IP(i),Mass.omega_IP(i)) = Wn;
        DAE.Fx(Mass.omega_IP(i),Mass.delta_HP(i)) = K12(i)*iMip(i);
        DAE.Fx(Mass.omega_IP(i),Mass.delta_IP(i)) = (-K12(i)-K23(i))*iMip(i);
        DAE.Fx(Mass.omega_IP(i),Mass.omega_IP(i)) = -(Dip(i)+D12(i)+D23(i))*iMip(i);
        DAE.Fx(Mass.omega_IP(i),Mass.delta_LP(i)) = K23(i)*iMip(i);
        DAE.Fx(Mass.delta_LP(i),Mass.omega_LP(i)) = Wn;
        DAE.Fx(Mass.omega_LP(i),Mass.delta_IP(i)) = K23(i)*iMlp(i);
        DAE.Fx(Mass.omega_LP(i),Mass.delta_LP(i)) = (-K23(i)-K34(i))*iMlp(i);
        DAE.Fx(Mass.omega_LP(i),Mass.omega_LP(i)) = -(Dlp(i)+D23(i)+D34(i))*iMlp(i);
        DAE.Fx(Mass.omega_LP(i),Syn.delta(i))     = K34(i)*iMlp(i);
        DAE.Fx(Syn.omega(i),Mass.delta_LP(i))     = K34(i)*iM(i);
        DAE.Fx(Syn.omega(i),Syn.delta(i))         = ...
            DAE.Fx(Syn.omega(k),Syn.delta(k)) + (-K34(i)-K45(i)).*iM(i);
        DAE.Fx(Syn.omega(i),Syn.omega(i))         = ...
            DAE.Fx(Syn.omega(k),Syn.omega(k)) - (D34(i)+D45(i)).*iM(i);
        DAE.Fx(Syn.omega(i),Mass.delta_EX(i))     = K45(i)*iM(i);
        DAE.Fx(Mass.delta_EX(i),Mass.omega_EX(i)) = Wn;
        DAE.Fx(Mass.omega_EX(i),Syn.delta(i))     = K45(i)*iMex(i);
        DAE.Fx(Mass.omega_EX(i),Mass.delta_EX(i)) = -K45(i)*iMex(i);
        DAE.Fx(Mass.omega_EX(i),Mass.omega_EX(i)) = -(Dex(i)+D45(i))*iMex(i);
        DAE.Fx(Mass.omega_HP(i),Mass.omega_IP(i)) = D12(i)*iMhp(i);
        DAE.Fx(Mass.omega_IP(i),Mass.omega_HP(i)) = D12(i)*iMip(i);
        DAE.Fx(Mass.omega_IP(i),Mass.omega_LP(i)) = D23(i)*iMip(i);
        DAE.Fx(Mass.omega_LP(i),Mass.omega_IP(i)) = D23(i)*iMlp(i);
        DAE.Fx(Mass.omega_LP(i),Syn.omega(k)) = D34(i)*iMlp(i);
        DAE.Fx(Syn.omega(k),Mass.omega_LP(i)) = D34(i)*iM(i);
        DAE.Fx(Syn.omega(k),Mass.omega_EX(i)) = D45(i)*iM(i);
        DAE.Fx(Mass.omega_EX(i),Syn.omega(k)) = D45(i)*iMex(i);
    end

end