应用介绍

电力系统的稳定性和控制，本程序的目的是分析小信号稳定性， Example 12.3 (page 752)，参考文档： http://www.apollocode.net/a/1535.html

  clear,clc,close all; % clear memory and workspace

s = tf('s'); % specifies the transfer function H(s) = s (Laplace variable)

digits(4);  % define precision, useful for symbolic math operations

deg2rad = pi/180;      % conversion factor from degs to radians

rad2deg = 1/deg2rad;   % conversion factor from radians to deg

%% Objectives:

% The objective of this example is to analyze the small-signal stability

% characterics of the system about the steady-state opeating condition

% following the loss of circuit 2.

%% Assumptions:

% The effects of amortisseurs may be neglected

% Neglect saliency of rotor

%% Operating Conditions:

Pt = 0.9;       %[w] active power supplied by generator

Qt = 0.3;       %[var] reactive power supplied by generator, overexcited)

Et_mag = 1;     %[Vll] magnitude, line-to-line terminal voltage

Et_angle = 36;  %[deg] 

Et=Et_mag*exp(j*deg2rad*(Et_angle));      % polar notation

% NOTE to get Et in rectangular form: abs(Et), rad2deg*(angle(Et))

EB_mag = 0.995;     %{Vll] magnitudeline-to-line bus voltage

EB_angle = 0;   %[deg]  

EB=EB_mag*exp(j*deg2rad*(EB_angle));      % polar notation

%% A 555MVA, 60Hz turbine generator has the following parameters:

% GIVEN SPECIFICAITON:

S = 555e6;      %[VA] total power base

f0 = 60;         %[Hz] frequency

Xd=1.81;    % note that in PU system, Ld = Xd

Xq=1.76;

Xp_d=0.3;   %X'd

Xl=0.16;

Ra = 0.003;

Tp_d0 = 8;  %[s]

H=3.5;

KD=0;

% The above parameters are unsaturated values

% The effect of saturation is to be represented by assuming that d and q

% axes have similar saturation characteristics with:

Asat =0.031; Bsat = 6.93; YTI = 0.8;  % based on definitions of section 3.8.2

% NOTE: giving the saturation characteristics in this manor means you have

% to CALCULATE saturation coefficients Ksd, Ksq.

%% Network Connection

% Generator is connected to an infinite bus through a step-up transformer

% (Xtrans = j0.15) and two parallel transmission lines (Xl1 = j0.5, Xl2 = 0.093)

% assuming lossless conditions.

Xtrans = 0.15;    %[ohm] impedance of step-up transformer

Xl1 = 0.5;        %[ohm] impedance of line 1 

Xl2 = 0.93;       %[ohm] impedance of line 2 

% Compute the equivalent network (single line)

Xe = Xtrans + Xl1;  %[ohm] equivalent network impedance, recall line 2 is out of service

Re = 0;             %[ohm] assume no losses

%% The per unit fundamental parameters

% elements of d- and q-axis equivalent circuits of the equivalent generator

% 

% First compute the unsaturated mutual inductances (pg. 154)

Ladu = Xd - Xl;     % mutual inductance, unsaturated value

Laqu = Xq - Xl;     % mutual inductance, unsaturated value

% Compute the rotor leakage inductances from the expressions for transient

% and subtransient inductances (pg. 154).  

% Based on equation 4.29, the expression for L'd based on the classical

% definition is:  L'd = Ll + Ladu*Lfd/(Ladu+Lfd); 

% note L'd =X'dm & Ll = Xl  -- L'd is a GIVEN VALUE (and is the unsaturated value)

% we need to solve for Lfd (unsaturated rotor field inductance)

% therefore:  L'd(Ladu+Lfd) = Ll*(Ladu + Lfd)+Ladu*Lfd

%             Lfd(L'd-Ll-Ladu) = -L'd*Ladu + Ll*Ladu

Lfd = (-Xp_d*Ladu + Xl*Ladu)/(Xp_d-Xl-Ladu);  % unsaturated value, rotor field inductance

Lp_adu = 1/((1/Ladu)+1/(Lfd)); % L'ad = X'ad mutual unsaturated value   %%%%%%%%  is this even used anywhere? %%%%%%%

%

% Next we compute the rotor resistances from the expressions for the

% open-circuit time constants: Equ. 4.15 and 4.23

% T'd0 = (Ladu+Lfd)/Rfd  [pu]

% NOTE: The time constant in per unit is equal to 377 times the time

% constant in seconds (pg 155)

Rfd = (Ladu + Lfd)/(377*Tp_d0);  %[pu]

%% Initial Steady-State Values:

St = Pt + j*Qt;     %[VA] total power at generator output in rect. notation

St_mag = abs(St);    %[VA] magnitude of total power 

St_angle = rad2deg*(angle(St));  %[deg] angle 

It = conj(St)/conj(Et);        % [A] phase current.  Take conjugate to recognize that current is positive out of generator

It_mag = abs(It);

It_angle = rad2deg*(angle(It));  % [deg] phase current angle  

% 

% SOLVE for saturated coefficients of Ksq, Ksd 

Yat0 = abs(Et + It*(Ra+j*Xl));

YIO = Asat*exp(Bsat*(Yat0-YTI));

Ksd = Yat0/(Yat0+YIO);  % saturation coefficient, same as value given on pg. 753

Ksq = Ksd;              % saturation coefficient, same as value given on pg. 753

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Find the saturated values for parameters

Lads = Ksd * Ladu;  % Ld mutual inductance, saturated value

Laqs = Ksq * Laqu;  % Lq mutual inductance, saturated value

Lp_ads = 1/((1/Lads)+(1/Lfd));  % saturated value inductance of L'd  (see eq. 12.91)

Lds = Lads+Xl;  % the SATURATED euiv. of Ld or Xd (previously defined)

Lqs = Laqs+Xl;  % the SATURATED euiv. of Lq or Xq (previously defined)

%%%Lp_ds = 1/((1/Lds)+(1/Lfd));   ???????????

Xads = Lads;    % equivalenet when in per unit

Xds = Lds;      % equivalenet when in per unit

Xqs = Lqs;      % equivalenet when in per unit

%Xp_ds = Lp_ds;   ???????????

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The network constraint equation 

% As a check, we know Ebus (EB) as it was given, 0.9950 <0.  

% Since, Et = EB + (RE + jXE)*It

EB_check = Et - (Re + j*Xe)*It;

EB_mag_check = abs(EB_check);

EB_angle_check = rad2deg*(angle(EB_check));

%

% Calculate the internal rotor angle (Si).  See Figure 3.23, pg 101.

%Check calculation of Si_angle using test data from text page 103

%Xqs = 1.494; It_mag = 1; Et_mag = 1; St_angle = 26 % test #'s page 103

Si = rad2deg*(atan((Xqs*It_mag*cos(deg2rad*(St_angle))-Ra*It_mag*sin(deg2rad*(St_angle)))/(Et_mag+Ra*It_mag*cos(deg2rad*(St_angle))+Xqs*It_mag*sin(deg2rad*(St_angle)))));

% Another way of calculating the internal rotor angle:

% Decompose into d-q components, see figure 12.6

Etq = Et + (Ra + j*Xqs)*It;

Etq_mag = abs(Etq);

Etq_angle = rad2deg*(angle(Etq));

Si_angle_verify = Etq_angle - Et_angle;  % [deg] internal angle (see Section 3.6.3), between Et and q-axis  

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Calculate Initial Conditions, as outlined on page 746, solutions given on pg. 753

ed0 = real(Et_mag*sin(deg2rad*(Si))); % real part of |Et|*sin(Si)

eq0 = real(Et_mag*cos(deg2rad*(Si))); % real part of |Et|*cos(Si)

id0 = real(It_mag*sin(deg2rad*(Si+St_angle))); % real part of |It|*sin(Si+phi)

iq0 = real(It_mag*cos(deg2rad*(Si+St_angle))); % real part of |It|*cos(Si+phi)

% 

% Solve for S0, equations on page 746, solutions given on pg. 753 

EBd0 = ed0 - Re*id0+Xe*iq0;

EBq0 = eq0 - Re*iq0-Xe*id0;

S0 = rad2deg*(atan(EBd0/EBq0));  

%

% Solve for Efd0, equations on page 746, solutions given on pg. 753 

%EB2 = sqrt((EBd0^2+EBq0^2));  % NOTE: this is same value as given for EB

ifd0 = (eq0+Ra*iq0+Lds*id0)/Lads;   % as per pg 746

Efd0 = Ladu*ifd0;                   % as per pg 746

Yad0 = Lads*(-id0+ifd0);             % as per pg 746

Yaq0 = -Laqs*iq0;                   % as per pg 746

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Small Signal Analysis

% "Since we are expressing small-signal performance in terms of perturbed

% values of flux linkages and currents, a distinction has to be made etween

% TOTAL saturation "s" and INCREMENTAL saturation "i".  Incremental

% saturation is associated with perturbed values of flux linkages and

% currents."

% Solve for Incremental values, Ksd(incr) & Ksq(incr)

Ksdi = 1/(1+Bsat*Asat*exp(Bsat*(Yat0 - YTI))); % incremental saturation, verify on page 753

Ksqi = Ksdi; 

Ladi = Ksdi*Ladu;   % Incremental mutual inductances

Laqi = Ksdi*Laqu;   % Incremental mutual inductances

Ldi = Ladi + Xl;    % Incremental inductance

Lqi = Laqi + Xl;    % Incremental inductance

Lp_adi = 1/((1/Ladi)+(1/Lfd));  % Incremental mutual inductance of L'd    

Xadi = Ladi;        % Equivalent in per unit system

Xaqi = Laqi;        % Equivalent in per unit system

Xdi = Ldi;          % Equivalent in per unit system

Xqi = Lqi;          % Equivalent in per unit system

Xp_di = Lp_adi;     % Equivalent in per unit system

%

%

% Preparing to linearize.. see page 742

% NOTE:  use of "incremental" values, as per pg.745

RT = Ra+Re;

XTq = Xe+Xqi;       % Text shows Xqs .. but here we use "i", the incremental value

XTd = Xe+Xp_di+Xl;  % Text shows Xqs .. but here we use "i", the incremental value

D = RT^2+XTq*XTd;

% Linearized system equations (pg. 742)

m1 = (EB_mag*(XTq*sin(deg2rad*(S0)) - RT*cos(deg2rad*(S0))))/D;

n1 = (EB_mag*(RT*sin(deg2rad*(S0)) + XTd*cos(deg2rad*(S0))))/D;

m2 = XTq/D * Ladi/(Ladi + Lfd);   

n2 = RT/D * Ladi/(Ladi+Lfd);      

% %%%% SUMMARY: Initial Steady-state Values of the System Varaibles

% disp('SUMMARY: Initial Steady-state Values of the System Varaibles:')

% Ksd

% Ksq

% Si

% ed0

% eq0

% id0

% iq0

% S0

% Efd0

% Ksdi

% Ksqi

% m1

% m2

% n1

% n2

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The Constants associated with the block diagram (fig. 12.9)

K1 = n1*(Yad0+Laqi*id0)-m1*(Yaq0+Lp_adi*iq0); % eq. 12.113  NOTE: incremental values used

K2 = n2*(Yad0 + Laqi*id0)-m2*(Yaq0+Lp_adi*iq0)+Lp_adi/Lfd*iq0; % eq. 12.114 NOTE: incremental values used

%% Construct State-Space Matrix

a11 = -KD/(2*H);

a12 = -K1/(2*H);

a13 = -K2/(2*H);

a21 = 2*pi*f0;

a22 = 0;

a23 = 0;

a31 = 0;

a32 = -((2*pi*f0)*Rfd/Lfd)*m1*Lp_adi;  % NOTE: use of incremental value

a33 = -((2*pi*f0)*Rfd/Lfd)*(1-Lp_adi/Lfd + m2*Lp_adi);  % NOTE: use of incremental value

% Construct the A-matrix

A = [a11 a12 a13; a21 a22 a22; a31 a32 a33]; % A-matrix of SMIB

%disp('The A-matrix:'); disp(A); %Displays the A-matrix to workspace

%

b11 = 1/(2*H);

b12 = 0;

b21 = 0;

b22 = 0;

b31 = 0;

b32 = (2*pi*f0)*Rfd/Ladu;   %%%%%  Why is this unsaturated mutual inductance Ld?

% The Constants associated with the block diagram (fig. 12.9)

K3 = -b32/a33;

K4 = -a32/b32;

T3 = -1/a33;

%% Eigenvalues 

% %%%%%%%%%%%%%%%%%%%%%%%%%  TEST from pg 734  %%%%%%%%%%%%%%%%%%%%%%%% 

A = [(-1.43) -0.108; 377 0]; % with KD  = 10  %%%%%%%%%%%%%%%%%%%%%%%%%%%

%

[eigen_Vec, eigen_Val] = eig(A);     % eigen vectors and eigenvalues, where the eigen values are in a diagonal matrix

disp('The Eignvalues:'); disp(diag(eigen_Val));  % Displays eigenvalues to workspace

%disp('Characteristic Equation:'); disp(poly(A));    % Displays characteristic equation

syms h; % define a symbol for lamda% At = [(-1.43) -0.108; 377 0];

Dh = diag(h)*eye(size(A));  % form a diagonal matrix of lambda, same size as A-matrix

dec= vpa((A-Dh));   % convert to decimal values

eq1 = sort(det(dec));     % characteristic equation is the determinant of (A-D), sorted 

disp('Characteristic Equation:')

pretty(eq1)   % displays result in a "pretty" way

%

% Define 2nd order characteristic equation variables

% syms phi wn;

% char_eq = h^2+2*phi*wn*h+wn^2;

% disp('of the form:')

% pretty(char_eq)

%

% Look at the characteristic equation enter in wn

% wn = sqrt(8.318);   % [rad] undamped natural frequency 

% fn = wn/(2*pi);     % [Hz] undamped natural frequency 

% sigma = 1.430/(2*wn);   % Damping ratio

% wd = wn*sqrt(1-sigma^2);    %[rad] damped frequency

% fd = wd/(2*pi);     % [Hz] damped frequency

% Find Right Eigenvector and eigenvalue, check with page 735

syms phi_11 phi_21 phi_31 phi_12 phi_22 phi_32 phi_31 phi_32 phi_33;  % manually determined

 [right,ev]=eig(A);

 disp('The Right Eigenvector Matrix:')

 disp(right)

% Method 2

% For the first eigenvalue

eig1= diag(eigen_Val(1))*eye(size(A));  % (lamda * I)

eig1a= vpa((A-eig1));   % (A-lamda*I) answer convert to decimal values

phi_m1 = [phi_11; phi_21;];  % first matrix of phi variables

eig1b = eig1a * phi_m1;

% one of the eigenvecotors corresesponding to an eigenvalue has to be set

% arbitrarily, therefore let phi_21 = 1, and solve for phi_11.

%phi_21 = 1;

[phi_11]=solve(subs(eig1b(1), phi_21, 1))

%

% Similarily, for the second eigenvalue

eig2= diag(eigen_Val(2))*eye(size(A));  % (lamda * I)

eig2a= vpa((A-eig2));   % (A-lamda*I) answer convert to decimal values

phi_m2 = [phi_12; phi_22;];  % first matrix of phi variables

eig2b = eig2a * phi_m2;

% one of the eigenvecotors corresesponding to an eigenvalue has to be set

% arbitrarily, therefore let phi_22 = 1, and solve for phi_12.

%phi_22 = 1;

[phi_12]=solve(subs(eig1b(2), phi_22, 1))

%

% Find left eigenvector, check with page 736

 left=inv(right);

 disp('The Left Eigenvector Matrix:')

 disp(left)

%% Plotting

% Plot Eigenvalues:

% figure(1)

% real_eigen=real(eig(A));

% imag_eigen=imag(eig(A));

% plot(real_eigen,imag_eigen,'*')

% xlabel ('Real')

% ylabel ('Imaginary')

% title ('System Eigen Values')

% grid

break

Identity = eye(size(A)); % returns an identity matrix the same size as A

Use [W,D] = eig(A.'); W = conj(W) %to compute the left eigenvectors.

参考书籍《电力系统的稳定性和控制》Prabha Kundur