Как составить матрицы состояний

Пространство состояний в задачах проектирования систем оптимального управления

Введение

Исследование системы управления во временной области с помощью переменных состояния широко используется в последнее время благодаря простоте проведения анализа.

Состоянию системы соответствует точка в определённом евклидовом пространстве, а поведение системы во времени характеризуется траекторией, описываемой этой точкой.

При этом математический аппарат включает готовые решения по аналоговому и дискретному LQR и DLQR контролерам, фильтра Калмана, и всё это с применением матриц и векторов, что и позволяет записывать уравнения системы управления в обобщённом виде, получая дополнительную информацию при их решении.

Целью данной публикации является рассмотрение решения задач проектирования систем оптимального управления методом описания пространства состояний с использованием программных средств Python.

Теория кратко

Векторно-матричная запись модели линейного динамического объекта с учетом уравнения измерения принимает вид:

(1)

Если матрицы A(t), B(t) и C(t) не зависят от времени, то объект называется объектом с постоянными коэффициентами, или стационарным объектом. В противном случае объект будет нестационарным.

При наличии погрешностей при измерении, выходные (регулируемые) сигналы задаются линеаризованным матричным уравнением:

(2)

где y(t) – вектор регулируемых (измеряемых) величин; C(t) – матрица связи вектора измерений с вектором состояний; v(t) – вектор ошибок измерений (помехи).

Структура линейной непрерывной системы, реализующая уравнения (1) и (2), приведена на рисунке:

Данная структура соответствует математической модели объекта, построенной в пространстве состояний его входных x(t), u(t), выходных y(t) и внутренних, или фазовых координат x(t).

Для примера рассмотрим математическую модель двигателя постоянного тока с независимым возбуждением от постоянных магнитов. Система уравнений электрической и механической частей двигателя для рассматриваемого случая будет выглядеть так:

(3)

Первое уравнение отражает взаимосвязь между переменными в цепи якоря, второе — условия механического равновесия. В качестве обобщенных координат выберем ток якоря I и частоту вращения якоря ω.

Управлением являются напряжение на якоре U, возмущением — момент сопротивления нагрузки Mc. Параметрами модели являются активное сопротивление и индуктивность цепи и якоря, обозначенные соответственно Rя, и Lя, а также приведенный момент инерции J и конструктивные постоянные се и см (в системе СИ: Cе=См).

Разрешая исходную систему относительно первых производных, получим уравнения двигателя в пространстве состояний.

(4)

В матричном виде уравнения (4) примут вид (1):

(5)

где вектор обобщенных координат , вектор управлений U =u (в рассматриваемом случае он является скаляром), вектор (скаляр) возмущений Mc=f. Матрицы модели:

(6)

Если в качестве регулируемой величины выбрать частоту вращения, то уравнение измерения запишется в виде:

(7)

а матрица измерений примет вид:

C=(0 1)

Сформируем модель двигателя в Python. Для этого вначале зададим конкретные значения параметров двигателя: U = 110 В; R =0,2 Ом; L = 0,006 Гн; J =0,1 кг/м2;Ce =Cm=1,3 В/С и найдем значения коэффициентом матриц объекта из (6).

Разработка программы формирующей модель двигателя с проверкой матриц на наблюдаемость и управляемость:

При разработке программы использовалась специальная функция def matrix_rank для определения ранга матрицы и функции, приведенные в таблице:

# -*- coding: utf8 -*-
from numpy import*
from control.matlab import *
from numpy.linalg import svd
from numpy import sum,where 
import matplotlib.pyplot as plt
def matrix_rank(a, tol=1e-8):
         s = svd(a, compute_uv=0)
         return sum( where( s>tol, 1, 0) )
u=110 # Напряжение якоря
J=.1  # Момент инерции
c=1.3; # Конструктивный коэффициент
R=.2; L=.006 # Активное сопротивление и индуктивность якоря
A=matrix([[-R/L ,-c/L],[ c/J ,0]])
print ('Матрица А : n %s'%A)
B=matrix([[1/L],[0]])
print ('Матрица B : n %s '%B)
C=matrix([[0, 1]])
print ('Матрица C : n %s '%C)
D=0
print ('Скаляр D : n %s '%D)
sd=ss(A,B,C,D) #Задание модели объекта в пространстве состояний
wd=tf(sd) # Задание передаточной функции двигателя
print ('Передаточная функция двигателя : n %s '%wd)
a=ctrb(A,B)
print(' Ранг матрицы управляемости : %s'%matrix_rank(a, tol=1e-8))
a=obsv(A,C)
print('Ранг матрицы наблюдаемости: %s'%matrix_rank(a, tol=1e-8))
y,x=step(wd) # Построение переходной характеристики  
plt.plot(x,y,"b")
plt.title('Переходная характеристика двигателя ') 
plt.ylabel('Частота вращения вала в рад/с')
plt.xlabel('Время С')
plt.grid(True)
plt.show()

Результаты работы программы:

Матрица А:
[[ -33.33333333 -216.66666667]
[ 13. 0. ]]
Матрица B:
[[166.66666667]
[ 0. ]]
Матрица C:
[[0 1]]
Скаляр D:
0
Передаточная функция двигателя:
2167/(s^2 + 33.33 s + 2817)
Ранг матрицы управляемости: 2
Ранг матрицы наблюдаемости: 2

Выводы:

1. На примере двигателя постоянного тока с независимым магнитным возбуждением рассмотрена методика проектирования управления в пространстве состояний;

2. В результате работы программы получены передаточная функция, переходная характеристика, а так же ранги матриц управляемости и наблюдаемости. Ранги совпадают с размерностями пространства состояний, что подтверждает управляемость и наблюдаемость модели.

Пример проектирования оптимальной системы управления с дискретным dlqr контролером и полной обратной связью

Определения и терминология

Линейно-квадратичный регулятор (англ. Linear quadratic regulator, LQR) — в теории управления один из видов оптимальных регуляторов, использующий квадратичный функционал качества.

Задача, в которой система описывается линейными дифференциальными уравнениями, а показатель качества, представляет собой квадратичный функционал, называется задачей линейно-квадратичного управления.

Широкое распространение получили линейно-квадратичные регуляторы (LQR) и линейно-квадратичные гауссовы регуляторы (LQG).

Приступая к практическому решению задачи всегда нужно помнить об ограничениях

Для синтеза оптимального дискретного регулятора линейных стационарных систем нужна функция численного решения уравнения Беллмана.Такой функции в библиотеке Python Control Systems [1] нет, но можно воспользоваться функцией для решения уравнения Риккати, приведенной в публикации [2]:

def dlqr(A,B,Q,R):
    """Solve the discrete time lqr controller.     
    x[k+1] = A x[k] + B u[k]     
    cost = sum x[k].T*Q*x[k] + u[k].T*R*u[k]
    """
    #ref Bertsekas, p.151 
    #first, try to solve the ricatti equation
    X = np.matrix(scipy.linalg.solve_discrete_are(A, B, Q, R))     
    #compute the LQR gain
    K = np.matrix(scipy.linalg.inv(B.T*X*B+R)*(B.T*X*A))     
    eigVals, eigVecs = scipy.linalg.eig(A-B*K)     
    return K, X, eigVals

Но нужно ещё учесть ограничения на синтез оптимального регулятора, приведенные в [3]:

• система, определяемая матрицами (A, B) должна быть стабилизируема;
• должны выполняться неравенства S> 0, Q – N/R–N.T>0, пара матриц (Q – N/R–N.T,
  A – B/R–B.T) не должна иметь наблюдаемые моды с собственными значениями на
  действительной оси.

После копаний в обширной и не однозначной теории, которую, по понятным причинам, я не привожу, задачу удалось решить, и все ответы можно прочитать прямо в комментариях к коду.

Структурная схема регулятора системы управления с обратной связью по всем переменным состояния изображена на рисунке:

Для каждого начального состояния x0 оптимальный линейный регулятор порождает оптимальное программное управление u*(x, k) и оптимальную траекторию х*(k).

Программа, формирующая модель оптимального управления с dlqr контролером

from numpy import *
import scipy.linalg
import matplotlib.pyplot as plt
def dlqr(A,B,Q,R):
          #Дискретное решение  матричного уравнения Риккати  x[k+1] = A x[k] + B u[k]
         P= matrix(scipy.linalg.solve_discrete_are(A, B, Q, R))
          #Дискретный контроллер DLQR
         K = matrix(scipy.linalg.inv(B.T*P*B+R)*(B.T*P*A))
         E, E1 = scipy.linalg.eig(A-B*K)
         return K, P, E
"""Параметры системы"""
A=matrix([[1, 0],[ -2 ,1]])
B=matrix([[1, 0],[1, 0]]).T
"""Параметры качества управления"""
Q=matrix([[0.5, 0],[0, 0.5]])
R=matrix([[0.5,0],[0, 0.5]])
T =100# Время регулирования
SS=0.5# Величина шаг
N =int(T/SS)# Количество шагов
K, P, E=dlqr(A,B,Q,R)#Вычисление параметров контроллера
print("K= n%s"%K)
print("P= n%s"%P)
print("E= n%s"%E)
x = zeros([2, N])
u= zeros([2, N-1])
x [0,0]=2
x [1,0]=1
for i in arange(0,N-2):
         u[:, i]= dot(- K,x[:, i])
         x[:, i+1]=dot(A,x[:, i])+dot(B,u[:, i])
x1= x[0,:] 
x2= x[1,:] 
t = arange(0,T,SS)
t1=arange(0.5,T,SS)
plt.subplot(221)          
plt.plot(t, x1, 'b')
plt.grid(True)
plt.subplot(222)
plt.plot(t, x2, 'g')
plt.grid(True)
plt.subplot(223)
plt.plot(t1, u[0,:], 'y')
plt.grid(True)
plt.subplot(224)
plt.plot(t1, u[1,:], 'r')
plt.grid(True)
plt.show()

Результаты расчета:

K=
[[ 0.82287566 -0.17712434]
[ 0.82287566 -0.17712434]]
P=
[[ 3.73431348 -1.41143783]
[-1.41143783 1.16143783]]
E=
[0.17712434+0.17712434j 0.17712434-0.17712434j]

Динамика состояний и управлений: x1, x2, u1, u2.

Вывод

Отдельные задачи оптимального управления по типу приведенных можно решать средствами Python, комбинируя возможности библиотек Python Control Systems, SciPy,NumPy, что, безусловно, способствует популяризации Python, учитывая, что ранее такие задачи можно было решать только в платных математических пакетах.

Ссылки

1. Python Control Systems Library.
2. Mark Wilfried Mueller LQR Controllers with Python.
3. Е.В.Никульчев. Пособие «Control System Toolbox»

ss

sys.InputName = 'controls';

sys.InputGroup.controls = [1 2];
sys.InputGroup.noise = [3 5];

sys.OutputName = 'measurements';

sys.OutputGroup.temperature = [1];
sys.OutputGroup.measurement = [3 5];

 sysarr.SamplingGrid = struct('time',0:10)

[zeta,w] = ndgrid(<6 values of zeta>,<9 values of w>)
M.SamplingGrid = struct('zeta',zeta,'w',w)

M(:,:,1,1) [zeta=0.3, w=5] =
 
        25
  --------------
  s^2 + 3 s + 25
 

M(:,:,2,1) [zeta=0.35, w=5] =
 
         25
  ----------------
  s^2 + 3.5 s + 25
 
...

A = [-1.5,-2;1,0];
B = [0.5;0];
C = [0,1];
D = 0;
sys = ss(A,B,C,D)

sys =
 
  A = 
         x1    x2
   x1  -1.5    -2
   x2     1     0
 
  B = 
        u1
   x1  0.5
   x2    0
 
  C = 
       x1  x2
   y1   0   1
 
  D = 
       u1
   y1   0
 
Continuous-time state-space model.

A = [0 1;-5 -2];
B = [0;3];
C = [0 1];
D = 0;

J = [8 -3 -3; -3 8 -3; -3 -3 8];
F = 0.2*eye(3);
A = -JF;
B = inv(J);
C = eye(3);
D = 0;
sys = ss(A,B,C,D)

sys =
 
  A = 
             x1        x2        x3
   x1  -0.04545  -0.02727  -0.02727
   x2  -0.02727  -0.04545  -0.02727
   x3  -0.02727  -0.02727  -0.04545
 
  B = 
           u1      u2      u3
   x1  0.2273  0.1364  0.1364
   x2  0.1364  0.2273  0.1364
   x3  0.1364  0.1364  0.2273
 
  C = 
       x1  x2  x3
   y1   1   0   0
   y2   0   1   0
   y3   0   0   1
 
  D = 
       u1  u2  u3
   y1   0   0   0
   y2   0   0   0
   y3   0   0   0
 
Continuous-time state-space model.

A = [-7,0;0,-10];
B = [5,0;0,2];
C = [1,-4;-4,0.5];
D = [0,-2;2,0];
ts = 0.2;
sys = ss(A,B,C,D,ts)

sys =
 
  A = 
        x1   x2
   x1   -7    0
   x2    0  -10
 
  B = 
       u1  u2
   x1   5   0
   x2   0   2
 
  C = 
        x1   x2
   y1    1   -4
   y2   -4  0.5
 
  D = 
       u1  u2
   y1   0  -2
   y2   2   0
 
Sample time: 0.2 seconds
Discrete-time state-space model.

A = [0 1;-5 -2];
B = [0;3];
C = [0 1];
D = 0;
Ts = 0.05;

sys = ss(A,B,C,D,Ts,'StateName',{'Position' 'Velocity'},...
    'InputName','Force');

A1 = [-1.5,-2;1,0];
B1 = [0.5;0];
C1 = [0,1];
D1 = 5;
sys1 = ss(A1,B1,C1,D1,'TimeUnit','minutes','InputUnit','minutes');

propValues1 = [sys1.TimeUnit,sys1.InputUnit]

propValues1 = 1x2 cell
    {'minutes'}    {'minutes'}

A2 = [7,-1;0,2];
B2 = [0.85;2];
C2 = [10,14];
D2 = 2;
sys2 = ss(A2,B2,C2,D2,sys1);

propValues2 = [sys2.TimeUnit,sys2.InputUnit]

propValues2 = 1x2 cell
    {'minutes'}    {'minutes'}

D = [2,4;3,5];
sys1 = ss(D)

sys1 =
 
  D = 
       u1  u2
   y1   2   4
   y2   3   5
 
Static gain.

H = [tf([1 1],[1 3 3 2]) ; tf([1 0 3],[1 1 1])];

State-space model with 2 outputs, 1 inputs, and 5 states.

sys = ss(H,'minimal');
size(sys)

State-space model with 2 outputs, 1 inputs, and 3 states.

load('identifiedModel.mat','ltiSys');

sysMeas = ss(ltiSys,'measured') 

sysMeas =
 
  A = 
            x1       x2
   x1    1.575  -0.6115
   x2        1        0
 
  B = 
        u1
   x1  0.5
   x2    0
 
  C = 
            x1       x2
   y1  -0.2851   0.3916
 
  D = 
       u1
   y1   0
 
  Input delays (sampling periods): 2 
 
Sample time: 0.04 seconds
Discrete-time state-space model.

sysNoise = ss(ltiSys,'noise')

sysNoise =
 
  A = 
           x1      x2      x3
   x1   1.026   -0.26  0.3899
   x2       1       0       0
   x3       0     0.5       0
 
  B = 
       v@y1
   x1  0.25
   x2     0
   x3     0
 
  C = 
             x1        x2        x3
   y1     0.319  -0.04738   0.07106
 
  D = 
          v@y1
   y1  0.04556
 
Input groups:        
    Name     Channels
    Noise       1    
                     
Sample time: 0.04 seconds
Discrete-time state-space model.

a = [2 -4; 4 2];
b = [-1; 0.5];
c = [-0.5, -2];
d = [-1];
e = [1 0; -3 0.5];
sysd = dss(a,b,c,d,e);

syse = ss(sysd,'explicit')

syse =
 
  A = 
        x1   x2
   x1    2   -4
   x2   20  -20
 
  B = 
       u1
   x1  -1
   x2  -5
 
  C = 
         x1    x2
   y1  -0.5    -2
 
  D = 
       u1
   y1  -1
 
Continuous-time state-space model.

bodeplot(sysd,syse,'g--')

a = realp('a',-1);
b = realp('b',3);

B = [-3.0;1.5];
C = [0.3 0];
D = 0;

Generalized continuous-time state-space model with 1 outputs, 1 inputs, 2 states, and the following blocks:
  a: Scalar parameter, 2 occurrences.
  b: Scalar parameter, 2 occurrences.

Type "ss(sys)" to see the current value and "sys.Blocks" to interact with the blocks.

A = [-1.5,-2;1,0];
B = [0.5;0];
C = [0,1];
D = 0;
sys = ss(A,B,C,D,'InputDelay',0.5,'OutputDelay',2.5)

sys =
 
  A = 
         x1    x2
   x1  -1.5    -2
   x2     1     0
 
  B = 
        u1
   x1  0.5
   x2    0
 
  C = 
       x1  x2
   y1   0   1
 
  D = 
       u1
   y1   0
 
  Input delays (seconds): 0.5 
  Output delays (seconds): 2.5 
 
Continuous-time state-space model.

                A: [2x2 double]
                B: [2x1 double]
                C: [0 1]
                D: 0
                E: []
           Scaled: 0
        StateName: {2x1 cell}
        StatePath: {2x1 cell}
        StateUnit: {2x1 cell}
    InternalDelay: [0x1 double]
       InputDelay: 0.5000
      OutputDelay: 2.5000
        InputName: {''}
        InputUnit: {''}
       InputGroup: [1x1 struct]
       OutputName: {''}
       OutputUnit: {''}
      OutputGroup: [1x1 struct]
            Notes: [0x1 string]
         UserData: []
             Name: ''
               Ts: 0
         TimeUnit: 'seconds'
     SamplingGrid: [1x1 struct]

A = [-1.2,-1.6,0;1,0,0;0,1,0];
B = [1;0;0];
C = [0,0.5,1.3];
D = 0;
sys = ss(A,B,C,D);

sysFeedback = feedback(sys,1);
P = pole(sysFeedback)

P = 3×1 complex

  -0.2305 + 1.3062i
  -0.2305 - 1.3062i
  -0.7389 + 0.0000i

A = [-0.5,-0.1;1,0];
B = [1;0];
C = [0,1];
D = 0;
sys = ss(A,B,C,D)

sys =
 
  A = 
         x1    x2
   x1  -0.5  -0.1
   x2     1     0
 
  B = 
       u1
   x1   1
   x2   0
 
  C = 
       x1  x2
   y1   0   1
 
  D = 
       u1
   y1   0
 
Continuous-time state-space model.

wc = 0.75;
C2 = pidtune(sys,'PID2',wc)

C2 =
 
                       1              
  u = Kp (b*r-y) + Ki --- (r-y) + Kd*s (c*r-y)
                       s              

  with Kp = 0.513, Ki = 0.0975, Kd = 0.577, b = 0.344, c = 0
 
Continuous-time 2-DOF PID controller in parallel form.

AG = [-3,0.4,0.3;-0.5,-2.8,-0.8;0.2,0.8,-3];
BG = [0.4,0,0.3,0.2,0;-0.2,-1,0.1,-0.9,-0.5;0.6,0.9,0.5,0.2,0];
CG = [0,-0.1,-1;0,-0.2,1.6;-0.7,1.5,1.2;-1.4,-0.2,0];
DG = [0,0,0,0,-1;0,0.4,-0.7,0,0.9;0,0.3,0,0,0;0.2,0,0,0,0];
sysG = ss(AG,BG,CG,DG)

sysG =
 
  A = 
         x1    x2    x3
   x1    -3   0.4   0.3
   x2  -0.5  -2.8  -0.8
   x3   0.2   0.8    -3
 
  B = 
         u1    u2    u3    u4    u5
   x1   0.4     0   0.3   0.2     0
   x2  -0.2    -1   0.1  -0.9  -0.5
   x3   0.6   0.9   0.5   0.2     0
 
  C = 
         x1    x2    x3
   y1     0  -0.1    -1
   y2     0  -0.2   1.6
   y3  -0.7   1.5   1.2
   y4  -1.4  -0.2     0
 
  D = 
         u1    u2    u3    u4    u5
   y1     0     0     0     0    -1
   y2     0   0.4  -0.7     0   0.9
   y3     0   0.3     0     0     0
   y4   0.2     0     0     0     0
 
Continuous-time state-space model.

AK = [-0.2,2.1,0.7;-2.2,-0.1,-2.2;-0.4,2.3,-0.2];
BK = [-0.1,-2.1,-0.3;-0.1,0,0.6;1,0,0.8];
CK = [-1,0,0;-0.4,-0.2,0.3];
DK = [0,0,0;0,0,-1.2];
sysK = ss(AK,BK,CK,DK)

sysK =
 
  A = 
         x1    x2    x3
   x1  -0.2   2.1   0.7
   x2  -2.2  -0.1  -2.2
   x3  -0.4   2.3  -0.2
 
  B = 
         u1    u2    u3
   x1  -0.1  -2.1  -0.3
   x2  -0.1     0   0.6
   x3     1     0   0.8
 
  C = 
         x1    x2    x3
   y1    -1     0     0
   y2  -0.4  -0.2   0.3
 
  D = 
         u1    u2    u3
   y1     0     0     0
   y2     0     0  -1.2
 
Continuous-time state-space model.

feedin = [4 2];
feedout = [1 3 4];
sys = feedback(sysG,sysK,feedin,feedout,-1)

sys =
 
  A = 
           x1      x2      x3      x4      x5      x6
   x1      -3     0.4     0.3     0.2       0       0
   x2    1.18   -2.56    -0.8    -1.3    -0.2     0.3
   x3  -1.312   0.584      -3    0.56    0.18   -0.27
   x4   2.948  -2.929   -2.42  -0.452   1.974   0.889
   x5   -0.84   -0.11     0.1    -2.2    -0.1    -2.2
   x6   -1.12   -0.26      -1    -0.4     2.3    -0.2
 
  B = 
            u1       u2       u3       u4       u5
   x1      0.4        0      0.3      0.2        0
   x2    -0.44       -1      0.1     -0.9     -0.5
   x3    0.816      0.9      0.5      0.2        0
   x4  -0.2112    -0.63        0        0      0.1
   x5     0.12        0        0        0      0.1
   x6     0.16        0        0        0       -1
 
  C = 
           x1      x2      x3      x4      x5      x6
   y1       0    -0.1      -1       0       0       0
   y2  -0.672  -0.296     1.6    0.16    0.08   -0.12
   y3  -1.204   1.428     1.2    0.12    0.06   -0.09
   y4    -1.4    -0.2       0       0       0       0
 
  D = 
          u1     u2     u3     u4     u5
   y1      0      0      0      0     -1
   y2  0.096    0.4   -0.7      0    0.9
   y3  0.072    0.3      0      0      0
   y4    0.2      0      0      0      0
 
Continuous-time state-space model.

State-space model with 4 outputs, 5 inputs, and 6 states.

Содержание:

Анализ электрических цепей методом пространства состояний:

Все переменные величины, характеризующие динамическую систему G (рис. 9.1) или имеющие определенное к ней отношение, целесообразно разделить на три группы: 1) входные переменные или входные воздействия

Величины    предполагаются функциями времени и обозначают соответственно значения величин  в момент времени t. Как известно, электрические процессы в любой схеме определяются параметрами и функциями источников а также независимыми начальными условиями переменных при — токами в индуктивных элементах и напряжениями на емкостных элементах  которые должны быть известны или рассчитаны. Через  и

выражаются искомые переменные состояния цепи или системы. Они же определяют энергетическое состояние схемы. Поэтому в качестве переменных состояния целесообразно выбирать токи  и напряжения  Источники, включенные в схему, можно назвать входными переменными или вектором входа искомые величины — выходными переменными или вектором выхода  Для цепи с независимыми токами и напряжениями должны быть заданы независимых начальных условий.

Множество всех значений, которые может принять вектор входа  в момент времени t образует пространство входа системы. Аналогично множество всех значений, которые может принять вектор выхода в момент времени образует пространство выхода системы, и множество всех значений, которые может принять вектор состояния в момент времени t образует пространство состояний системы. В любой момент времени t состояние системы определяется вектором  являющимся функцией вектора начального состояния и вектора входа  т. е.

Вектор выхода также является функцией и может быть записан как

Уравнения (9.1) и (9.2) часто называют уравнениями состояния системы. Для систем, описываемых линейными дифференциальными уравнениями, уравнения состояния (9.1) и (9.2) сводятся к следующим:

где — матрица коэффициентов размерности — матрица управления (выхода) размерности — матрица выхода размерности — матрица обхода системы раз-мерности

На рис. 9.2 изображена обобщенная схема, динамика которой описывается уравнениями (9.3) и (9.4) [4, 8].

Если система стационарная, то ее динамика описывается уравнениями состояния, матрицы которых имеют элементы, не изменяющиеся во времени,  т. е.

Матрица коэффициентов А определяет структуру системы, параметры элементов и их взаимные связи. Динамические свойства системы в основном определяются этой матрицей. Матрица управления В показывает связь управляющих (возмущающих) воздействий переменными состояния. При отсутствии входных воздействий, когда требуется исследовать свободное движение системы, матрица управление В равна нулю. Матрица С формирует выходные переменные из переменных состояний И наконец, матрица обхода системы D показывает непосредственное влияние входных воздействий на выходные координаты без преобразования вектора динамической системой.

Следует отметить, что уравнения состояния (9.5) и (9.6) — это матричная запись системы линейных дифференциальных уравнений 1-го порядка с постоянными коэффициентами, описывающими динамику соответствующей линейной системы а взаимные связи между переменными состояния, входными и выходными величинами.

Уравнения (9.5) и (9.6) принято называть стандартной формой записи уравнений динамики линейных управляемых систем с постоянными параметрами, имеющих произвольную структуру и произвольное число входов и выходов. Стандартная форма отличается компактностью и удобством преобразования. Процедура решения уравнений состояния в конечном итоге сводится к матричным преобразованиям над А, В, С и D, что весьма удобно для программирования на цифровых ЭВМ. Уравнения (9.5) и (9.6) являются исходной информацией при исследовании и проектировании систем управления методом пространства состояний. К такому виду можно привести формы записи уравнений динамики (передаточные функции, дифференциальные уравнения высокого порядка, матричные передаточные функции многомерных систем и т. д.), применяемые в классических методах исследования систем.

Получение уравнений состояния является начальным этапом исследования систем и цепей в современной теории управления.

Методы составления уравнений состояния электрических цепей и динамических систем

Чтобы найти уравнения состояния (9.5) и (9.6), необходимо динамику электрических цепей представить системой дифференциальных уравнений 1-го порядка. В качестве иллюстрации рассмотрим систему третьего порядка, описываемую уравнением

Для записи этого уравнения в векторно-математической форме положим

и тогда вместо уравнения (9.7) с учетом (9.8) получим: 

Система уравнений (9.9) совместно с уравнением (9.10) в матричной форме запишется:

Уравнения (9.11) и (9.12) дают конкретный вид матриц А, В, С и D уравнений (9.5) и (9,6).

Приведенный выше метод может с успехом применяться для систем с одним входом и выходом. В многоконтурных системах с несколькими входами и выходами указанная процедура реализуется не так просто поэтому существуют другие способы получения уравнений состояния (9.5) и (9.6).

Наиболее распространенным способом получения уравнений состояния исследуемой системы является представление ее в виде схемы системы в переменных состояния. Эта схема составляется из интеграторов, усилителей и суммирующих устройств. Обычно выходы интеграторов выбираются в качестве координат (переменных) состояния системы. Схема в переменных состояния даст наглядную физическую интерпретацию координат системы и описывает их взаимную связь. Схемы непрерывных систем в переменных состояния совпадают со схемами моделирования этих систем на аналоговых вычислительных машинах. Существует много разновидностей схем моделирования для одной и той же системы [8], отсюда и неоднозначность описания системы управления уравнениями состояния.

Сначала рассмотрим методы построения схем в переменных состояния для систем с одним входом и выходом, динамика которых описывается дифференциальным уравнением вида

при начальных условиях  

Затем будут показаны некоторые примеры, применяемые при построении схем в переменных состояния для систем более сложной структуры, динамика которых не может быть описана одним уравнением вида (9.13). Конечным результатом построения схем в переменных состояния будет получение уравнений состояния (9.5) и (9.6).

Метод комбинирования производных

Рассмотрим сначала более простой случай, когда начальные условия уравнения (9.13) нулевые. При этом уравнение (9.13) в операторной форме запишется как

где  

Представим уравнение (9.14) в виде отношения 

Это уравнение заменим следующими двумя:

где — изображение вспомогательной переменной

В дифференциальной форме последние два уравнения примут вид

Таким образом, вместо решения уравнения (9.13) можно решить уравнения (9.15) и (9.16). Составим схему решения уравнения (9.15), для этого разрешим его относительно старшей производной:

Предположим, что -я производная функции известна, тогда, подавая ее на цепочку последовательно соединенных интеграторов, на выходе каждого последующего шифратора будем иметь производные от функции более низких порядков. Если каждую из этих производных умножить на соответствующий коэффициент а, и образовать их отрицательную сумму вместе с функцией  то получим правую часть уравнения (9.17), т. е. -ю производную функции которая ранее предполагалась известной. Функция  которая является решением уравнения (9.13), получается линейной комбинацией производных от взятых с коэффициентами в соответствии с уравнением (9.16). Схема решения уравнений (9.15) и (9.16) приведена на рис. 9.3.

При ненулевых начальных условиях уравнения (9.13) необходимо их перевести в начальные условия (9.15). Предполагая, что уравнения (9.15) и (9.16) эквивалентны (9.13) и при ненулевых начальных условиях, запишем эти уравнения в операторной форме:

где  — полиномы начальных условий.

Подставляя из последних двух уравнений (9.18) в первое, получаем

Из этого уравнения находится тождественное равенство, по которому определяются начальные условия уравнения (9.15) по известным начальным условиям (9.13)

Приравнивая коэффициенты при одинаковых степенях р этого равенства, находим начальные условия интеграторов схемы на рис. 9.3.

Если в качестве переменных состояния выбрать выходные величины интеграторов схемы (см. рис. 9.3), то получим следующую систему уравнений:

Записывая эту систему уравнений в матричной форме вида (9.5) и (9.6)

получаем конкретный вид матриц А, В, С и D.

Применение этого метода для составления уравнений состояния не требует преобразования уравнения (9.13). Можно непосредственно по виду уравнения (9.13) составлять схему в переменных состояния, так как его коэффициенты являются и коэффициентами схемы в переменных состояния. Пересчет начальных условий происходит по уравнению (9.19).

Пример 9.1.

Составить уравнения состояния для системы, динамика которой описывается дифференциальным уравнением

Схема в переменных состояния приведена на рис. 9.4. Из схемы в переменных состояния получаем уравнения:

Эта система уравнений в матричной форме имеет вид:

Следовательно, 

Метод последовательного интегрирования

Запишем уравнение (9.13) и операторной форме при нулевых начальных условиях в виде

откуда

Составим цепочку из последовательно соединенных интеграторов. Сигнал на выходе крайнего правого интегратора примем за . Выполнение равенства (9.22) обеспечится в том случае, соли на вход крайнего левого и интегратора подать разность на вход следующего — разность  и выход предыдущего интегратора и т. д. Схема решения уравнения (9.22) тогда примет вид, показанный на рис. 9.5.

Если снова выбрать в качестве переменных состояния выходы интеграторов схемы рис. 9.5, то получим систему уравнений:

При ненулевых начальных условиях уравнения (9.13) начальные условия интеграторов схемы определяются соотношением [8]

Записывая эту систему уравнений в форме уравнений (9.5) и (9.6), находим матрицы А, В, С и D:

Как и для метода комбинирования производных, коэффициенты уравнения (9.13) являются одновременно и коэффициентами схемы и переменных состояния. Поэтому, зная общую структуру схемы в переменных состояния, можно непосредственно по виду уравнения (9.13) построить соответствующую схему в переменных состояния и найти уравнения состояния в виде (9.5) и (9.6).

Пример 9.2. 

Составить уравнения состояния для системы, описываемой уравнением

Схема а переменных состояния, построенная методом непосредственно интегрирования, показана на рис. 9.6. Из схемы получаем систему уравнений 1-го порядка:

В матричной форме эта система уравнений имеет вид:

Из последних уравнений ясно видны матрицы А, В, С. D.

В примере 9.2 матрица D отлична от нуля. Это имеет место в тех случаях, когда  т.е. максимальная степень производной правой части дифференциального уравнения равна максимальной степени производной левой части уравнения. Если то матрица D всегда будет нулевой матрицей. Случай соответствует физически нереализуемой системе, поэтому представляет только чисто теоретический интерес и здесь не рассматривается.

Метод разложения передаточной функции на элементарные дроби

Весьма перспективным для анализа цепей и систем методом пространства состояний является построение схемы в переменных состояния путем разложения передаточной функции на элементарные дроби. Суть его заключается в следующем. Уравнение (9.13) представляется в виде передаточной функции

Разложим эту передаточную функцию на элементарные дроби

Отсюда

Коэффициент d будет отличен от нуля при Схема в переменных состояния уравнения (9.26) показана на рис. 9.7. Уравнения состояния имеют вид:

Как видно из выражения (9.27), матрица А является диагональной.

Пример 9.3.

Составим уравнения состояния для электрической цепи второго порядка, передаточная функция которой имеет вид

Полюса этой передаточной функции равны:

Представим передаточную функцию в виде суммы элементарных дробей

Коэффициенты

Таким образом, необходимо составить схему в переменных состояния для уравнения

Схема в переменных состояния показана на рис. 9.8. Уравнения, связывающие переменные состояния  входную и выходную величины имеют вид

В матричном виде эти уравнения запишем таким образом:

При наличии кратных полюсов передаточной функции (9.25) матрица А будет представлена в канонической форме Жордана [8]. Ниже будет показано, что вычислительная процедура значительно упрощается, если матрица А диагональная или имеет каноническую форму Жордана. Основная трудность получения уравнений состояния вида (9,27) и (9,28) состоит в нахождении полюсов передаточной функции (9.25).

Преобразование неоднородных уравнений состояния в однородные

В результате построения схемы и переменных состояния можно получить уравнения состояния в виде (9.5). Это система неоднородных дифференциальных уравнений первого порядка. Их решение будет содержать две составляющие — свободную и вынужденную. Первая зависит от динамики системы и начального значения вектора переменных состояния  Вторая зависит также от матрицы А и вектора входных воздействий

Свободная составляющая определяется решением однородного уравнения

когда вектор входного воздействия  в выражении (9.5) предполагается нулевым.

Имеет смысл попытаться преобразовать систему неоднородных уравнений (9.5) в однородную вида (9.29). Цель этого преобразования заключается в том, чтобы получить уравнения состояния системы при наличии внешних воздействий , содержащие только переменные состояния. При этом значительно упрощается решение уравнений состояния.

Преобразование системы неоднородных уравнений в однородную можно осуществить в том случае, если вектор внешних воздействий представляется в виде решения определенного линейного однородного дифференциального уравнения (системы уравнений) с постоянными параметрами. Схема в переменных состояния для формирования вектора будет иметь свои переменные состояния с начальным значением

Рассмотрим методику построения схем в переменных состояния для формирования некоторых часто встречающихся входных воздействий

Построение схем в переменных состояния для типовых входных воздействий

Входное воздействие в виде полиномиальной функции:

Полиномиальное входное воздействие описывается функцией

при Обозначим

и продифференцируем уравнение (9.31) по t. Тогда

Аналогично

Продолжая этот процесс до тех пор, когда очередная производная будет равно нулю, получим следующую систему уравнений:

Величины называются переменными (координатами) состояния входа. Система дифференциальных уравнений (9.32) показывает, что полиномиальное воздействие можно получить с помощью цепочки последовательно соединенных интеграторов, как изображено на рис. 9.9.

Начальные условия на интеграторах:

Выходом этой системы является заданное входное воздействие, определяемое переменной состояния Матрицы уравнений состояния (9.32) находятся из выражений (9.32) и (9.31)

Пример 9.4. 

Составить схему в переменных состояния для формирования линейного воздействия и найти соответствующие уравнения состояния.

Систему уравнений, решением которой будет линейная функция, получим из выражения (9.32)

При начальных условиях выходной сигнал

Схема о переменных состояния показана на рис. 9.10,

В матричном виде уравнения состояния схемы (рис. 9.10) будут иметь вид
      

отсюда матрицы равны:

Входное воздействие в виде гармонической функции

В общем случае гармоническое воздействие представляется в виде

Решение дифференциального уравнения

даст гармонический сигнал. Это уравнение может быть записано в виде двух уравнений первого порядка

с начальными условиями

Схема в переменных состояния, которая следует из формул (9.37), показана на рис. 9.11. Сигнал определяется выражением

Изменяя начальные условия можно получать гармонические сигналы различной фазы и амплитуды.

Входное воздействие в виде экспоненциальной функции

Экспоненциальное входное воздействие описывают функцией

Пусть Тогда

Начальное значение Схема в переменных состояния системы, использующая это последнее уравнение, показана на рис. 9.12. Переменная является воздействием экспоненциального вида.

Методика преобразования неоднородных уравнений состояния в однородные

Для преобразования системы неоднородных уравнений (9.5) в однородную (9.29) необходимо вектор  рассматривать совместно с переменными состояния системы , т. е. включить вектор в число координат системы (рис. 9.13).

Допустим, что входной вектор есть решение уравнений состояния

полученных путем построения схемы и переменных состояния дифференциального уравнения, решением которого и является вектор Вектор будет входным для исследуемой системы с уравнениями состояния вида (9.5) и (9.6). Подставим уравнение (9.42) в (9.5) и (9.6)

Объединим в одну систему уравнения (9.43) и (9.41)

Обозначим расширенный вектор, включающий как переменные состояния управляемой системы, так и переменные состояния схемы в переменных состояния входного воздействии через т. е.

Уравнения (9.43), (9.44) и (9.45) объединим и запишем в виде:

Обозначая

получаем

Матрица побила название расширенной матрицы коэффициентов, — расширенной матрицы выхода, a — расширенного вектора состояния.

Применение методики преобразования неоднородной системы уравнений состояния в однородную проиллюстрируем на примере.

Пример 9.5.

Найти матрицы для системы примера 9.1 при линейном входном воздействии. Из примера 9.1

Для линейного воздействия

Подставляя в (9.49), получаем

На практике часто пользуются следующим приемом при преобразовании неоднородных уравнений (9,5) и (9.6) в однородные. Составляется расширенная схема в переменных состояния, включающая схемы в переменных состояния входного воздействия и собственно динамической системы. Из этой схемы непосредственно получают расширенные матрицы

Пример 9.6.

Для условий примера 9.5 найти матрицы (рис. 9.14).

Расширенная схема о переменных состояния показана на рис. 9.14. Из нее получаем:

В матричной форме эту систему уравнений можно записать, если обозначить  так

Результаты примеров 9.5 и 9.6 совпадают.

Формы решения уравнений состояния

Рассмотренные методы построения схем в переменных состояния позволяют получить математическое описание динамики линейной стационарной системы и электрической цепи в виде векторно-матричной системы дифференциальных уравнений первого порядка. Рассмотрим решение этих уравнений.

Форма решения однородных уравнений состояния

Если на систему не подаются внешние воздействия, то и уравнение (9.5) становится однородным

Решение этого уравнения описывает динамику системы за счет ненулевых начальных условий (свободное движение), когда внешние силы равны нулю. Предположим, что движение начинается в момент  из начального состояния Решение уравнения (9.52) можно получить по аналогии с решением скалярного уравнения в виде

Действительно, если подставить решение (9,53) в (9.52), предварительно взяв производную от (9.53), то получим тождество. Следовательно, (9.53) является решением однородного матричного уравнения (9.52).

Если обозначить

то уравнение (9.53) можно записать как

Выходной вектор системы будет иметь вид

Матрица называется матрицей перехода (переходной матрицей) состояния системы. Матрица перехода описывает движение конца вектора состояния в пространстве состояния из начальною положения а следовательно, и изменение (переход) состояния системы. Вектор дает возможность контролировать все координаты (переменные) состояния системы, поэтому решение (9.55) несет в себе полную информацию о динамике системы. Необходимо отметить, что объем вычислений при определении матрицы перехода состояния системы обычно больше, чем при решении линейного дифференциального уравнения относительно зависимой переменной. Однако имеющаяся дополнительная информация в результате вычисления матрицы перехода позволяет провести более глубокий и всесторонний анализ динамики исследуемой системы и электрических цепей.

Форма решения неоднородных уравнений состояния

Решение уравнения (9.5) будем искать в форме, аналогичной выражению (9.55). Положим

где — вектор, зависящий от времени и заменяющий вектор начального состояния в уравнении (9.55) при отсутствии внешних воздействий.

Дифференцируя выражение (9.57) по t, получаем

Учитывая, что уравнение (9.58) запишем в виде

Если выражение (9.57) является решением уравнения (9.5), то величины в правых частях уравнений (9.5) и (9.59) должны быть одинаковыми.

Отсюда

Решая это уравнение относительно получаем

Подставив выражения (9.61) в (9.57), имеем

Учитывая, что

выражение (9.62) примет вид

Постоянную интегрирования  можно определить, если в уравнении (9.64) принять Тогда

Так как то

Таким образом,

Вектор выхода будет иметь вид

Первое слагаемое выражения (9.67) — составляющая выходного вектора за счет ненулевых начальных условий  Она совпадает с решением однородного уравнения (9.52). Второе слагаемое характеризует реакцию системы на внешнее воздействие  третье — представляет собой непосредственное влияние вектора входа на выход

Аналитический подход к вычислению матрицы перехода

Формы решения однородной (9.56) и неоднородной (9.67) систем уравнений состояния содержат матрицу перехода Поэтому одной из центральных задач исследования систем является нахождение матрицы перехода. Существует много методов определения Основные подходы базируются на теореме Сильвестра [8, 33], а также к ним относится метод разложения  в бесконечный ряд и метод комплексной плоскости, которые мы и рассмотрим ниже.

Далее рассматриваются и обосновываются аналитические методы вычисления матрицы перехода применяемые в практике анализа и синтеза электрических цепей и динамических систем. Они иллюстрируются примерами.

Метод разложения матрицы перехода в бесконечный ряд

Переходная матрица раскладывается в бесконечный ряд

Этот метод наиболее трудоемок, если элементы матрицы «не исчезают» при небольших После выполнения суммирования следует найти в замкнутом виде каждый элемент матрицы  Это непростая задача, и если порядок недостаточно низкий, то задача может оказаться и непреодолимой. Однако метод разложения и бесконечный ряд весьма удобен при вычислении на цифровых вычислительных машинах, что и будет показано несколько ниже.

Пример 9.7.

Найти матрицу перехода для матрицы

Степени находят последовательным умножением

где  определяют по рекуррентной формуле

Переходная матрица

Если свернуть бесконечные ряды внутри матрицы получим

Метод комплексной плоскости

Возьмем преобразование Лапласа от обеих частей уравнения (9.52)

отсюда

Применив к уравнению (9.70) обратное преобразование Лапласа,  получим

Сравнивая выражения (9.53) и (9.71), приходим к выводу, что

Основная трудность этого метода состоит в нахождении матрицы, обратной

Пример 9.8.

Найти матрицу перехода , если

Найдем матрицу

Матрица, обратная имеет вид

Пользуясь таблицей преобразования Лапласа [26, 28], находим

Рассмотренный метод комплексной плоскости дает возможность проследить физическим смысл матрицы перехода Произвольный -й член уравнения (9.55) при можно записать в виде суммы

где — элемент матрицы стоящий на пересечении -й строки и -го столбца. Элемент  можно определить, положив равным единице и считая начальные условия для остальных переменных состояния равными нулю, тогда Применительно к схемам в переменных состояния это эквивалентно подаче единичного начального условия на выход интегратора и наблюдению за выходом интегратора  Однако единичное начальное условие на выходе интегратора равносильно подаче единичного мгновенного импульса на вход этого интегратора. Следовательно  представляет собой реакцию на выходе интегратора при подаче единичного мгновенного импульса на вход интегратора в то время как начальные условия остальных интеграторов равны нулю. Поэтому можно интерпретировать как импульсную переходную функцию от входа сумматора на входе интегратора к выходу интегратора Эти импульсные переходные функции образуют

При рассмотрении физического смысла матрицы перехода было указано, что каждый элемент можно интерпретировать как импульсную переходную функцию определенной части схемы в переменных состояния. Выражение (9.73) показывает, что реакция системы на единичный мгновенный импульс будет содержать конечное число составляющих которые получили название «частот» системы [8, 33]. Они позволяют проводить оценку характера решения уравнений состояния и дают возможность творчески решать задачи расчета и проектирования электрических цепей и систем. Например, путем выбора начальных условий или специальным формированием выходной координаты можно добиться отсутствия на выходе нежелательной «частоты» системы.

Построение алгоритмов вычисления матрицы перехода

Основой для разработки алгоритмов вычисления матрицы перехода   служат рассмотренные выше методы, применяемые для аналитического подхода к ее определению. Рассмотрим алгоритм по методу разложения в бесконечный ряд.

При использовании этого метода для машинного определения матрицы возникают две самостоятельные задачи. Во-первых, вычисление каждого последующего члена бесконечного ряда

с одновременным их суммированием и, во-вторых, оценка точности вычисления, поскольку ряд (9.74) необходимо ограничивать конечным числом членов.

Алгоритм вычисления бесконечного ряда

При решении этой задачи воспользуемся рекуррентными соотношениями, позволяющими наиболее просто и экономно составить алгоритм вычисления бесконечною ряда (9.74). Обозначим

Тогда

при

Таким образом, вычисление ряда (9.74) состоит из последовательных циклов, в каждом из которых будут использоваться результаты предыдущих вычислений. На рис. 9.15 показана схема алгоритма вычисления использующего рекуррентные соотношения (9.77) и (9.78). Алгоритм не имеет условия, по которому будет закончен счет ( ограничено), и поэтому число членов ряда (9.74) не определено.

 

Критерий ограничения бесконечного ряда

Вычисление матрицы перехода с точностью, не хуже заданной, требует способа оценки остаточного члена бесконечного ряда (9.74). Если ряд (9.74) ограничить первыми членами, то остаточный член примет вид

Представим остаточный член (9.79) в следующем виде:

Заменим правую часть ряда (9.80) мажорирующим матричным рядом, где каждый член будет не меньше соответствующего члена ряда (9.80). Тогда вместо (9.80) можно записать

Правая часть этого выражения является матричной геометрической прогрессией с начальным членом и знаменателем  Если норма знаменателя геометрической прогрессии меньше единицы, т. е.
(9.82)

то выражение (9.81) запишется так [8]:

Таким образом, остаточный член матричного ряда (9.74) можно оценивать по матрице в правой части выражения (9.83) при выполнении условия (9.К2). С увеличением каждый член матрицы (9.83) будет уменьшаться. Значение при котором наибольший член матрицы (9.83) будет меньше заданной точности и определяет число членов ряда (9.74).

Использование критерия (9.83) для ограничения числа членов бесконечного ряда требует вычисления нормы матрицы, так как

выражение (9.83) справедливо только при выполнении условия (9.82). На рис. 9.16 показана схема алгоритма вычисления нормы произвольной матрицы А. В качестве ее нормы выбрана максимальная сумма абсолютных значений элементов строки матрицы А. Этот алгоритм должен быть составной частью оценки числа членов ряда (9.74).

На рис. 9.17 приведена схема алгоритма, но которому можно определить число членов ряда (9.74), обеспечивающих точность вычисления матрицы перехода не хуже заданной. Алгоритм использует правую часть выражения (9.83) и имеет итерационный характер.

Алгоритм (рис. 9.17) требует достаточно сложной операции обращения матрицы. Целесообразно матричную оценку точности вычисления (9.83) заменить скалярной. В качестве такой оценки возьмем норму остаточного числа

Проведем преобразования правой части критерия (9.84)

Известно [30], что если

  и

  то

Заменяя норму обратной матрицы в выражении (9.85) и учитывая, что

получаем      

Рассмотрим построение алгоритма вычисления матрицы перехода с использованием скалярной оценки (9.90) остаточного члена ряда (9.74).

Алгоритм должен содержать цикл пои при каждом новом значении вычислять очередной член ряда (9.74) и норму матрицы Если

то можно использовать скалярную оценку остаточного члена матричного ряда (9.74). Для этого необходимо определять норму последующего слагаемого ряда (9.74) и при выполнении условия

ограничивать число членов ряда (9.74). Поскольку

всегда найдется такое значение которое обеспечит требуемую точность На рис. 9.18 приведена схема алгоритма вычисления матрицы перехода включающая скалярную оценку остаточного члена.

Ниже представлен вариант программы вычисления переходной матрицы  в соответствии с алгоритмом (рис. 9.18), составленный в среде Mathcad [43].

Исходные параметры системы (из примера 9.9.5):

Программа вычисления матрицы перехода

Метод разложения   в бесконечный ряд позволяет определить матрицу перехода для одного момента времени t. Если требуется иметь зависимость матрицы перехода   от времени, то необходимо вычислять для дискретных моментов времени, в результате этого получим функцию перехода в табличном виде. Такая информация может быть использована при анализе характера поведения во времени, а также при использовании решений уравнений состояния, представленных выражениями (9.53) и (9.62).

Иногда, особенно при решении задач проектирования, требуется иметь функцию перехода в аналитическом виде, т. е. в виде матрицы, элементами которой являются функции времени. Такой вид дает возможность более глубоко и с меньшими затратами времени оценить характер поведения и тем самым наметить пути совершенствования системы и повышения ее качественных показателей. Рассматриваемые ниже методы вычисления матрицы перехода позволяют решать эти задачи.

Аналитический подход к решению уравнений состояния

Формы решения уравнений состояния получены в разд. 9.4. Форма решения (9.53) применяется при исследовании как свободного движения системы, так и вынужденного, если входной вектор можно представить как решение некоторой линейной системы дифференциальных уравнений первого порядка, т.е. применить процедуру расширения матрицы коэффициентов и матрицы выхода. Форма решения (9.62) является более общей и, естественно, ее практическое использование требует больших усилий и навыков исследования. Приведенные ниже примеры показывают процедуру аналитического решения уравнений состояния: при исследовании свободного движения системы (пример 9.9), свободного и вынужденного движения при построении расширенной схемы в переменных состояния (пример 9.10).

Пример 9.9.

Найти свободное движение системы, схема которой в переменных состояния приведена на рис. 9.4, при начальных условиях Матрицы А, В, С и D для этой системы получены в примере 9.1

Поскольку требуется исследовать свободное движение, то а уравнения состояния становятся однородными, и их решение имеет вид (9.53)

где  из примера 9.7

Таким образом

Пример 9.10. Найти реакцию  системы 

и ошибку при линейном воздействии и начальных условиях

Расширенная схема в переменных состояния показана на рис. 9.19. Уравнения состояния расширенной системы имеют вид

при начальных условиях

Обозначая получаем

Решение этой системы уравнений можно записать в виде

.

Матрицу перехода представим бесконечным рядом

Степени найдем последовательным умножением:

Для нахождения вообще говоря, нет необходимости определять все элементы матрицы поскольку и вектор имеют нулевые элементы. Поэтому определим, какие из элементов матрицы необходимы при решении поставленной задачи. Представим в общем виде и найдем произведение

Для определения  достаточно знать два элемента матрицы 

Таким образом,

Входное воздействие связано с переменной состояния соотношением (см. рис. 9.19)

Переменная состояния

Учитывая, что получаем

Этот результат очевиден из условия задачи. Ошибка равна

На рис. 9.20 показаны При установившаяся ошибка равна 1,5.

Построение алгоритмов решения уравнений состояния

Рассмотренные примеры показывают основные вычислительные трудности, возникающие при решении задач исследования линейных динамических систем. Более общая форма решения (9.62) требует аналитического интегрирования достаточно сложных функций, особенно если вектор входных воздействий не может быть представлен комбинацией стандартных функций. Может возникнуть такая практическая ситуация, при которой вообще невозможно будет провести интегрирование. Кроме того, аналитический метод решения может применяться только при аналитически заданном векторе входных воздействий. Следовательно, решать задачи при табличном или графическом представлении невозможно. Все вышесказанное требует разработки алгоритмов решения уравнений состояния и использования вычислительных машин для исследования динамических систем. Сначала рассмотрим построение алгоритма для исследования свободного движения системы, который может без изменений применяться и для исследования вынужденного движения при наличии расширенных уравнений состояния. Затем будет описан общий случай, соответствующий алгоритмизации формы решения (9.62).

Построение алгоритмов решения однородных уравнений состояния

Простой алгоритм. Выходной вектор при отсутствии входных воздействий определяется по формуле (9.53)

Построение алгоритма для вычисления достаточно подробно было рассмотрено в разд. 9.6, поэтому в дальнейшем будем считать, что алгоритм определения известен, и на схеме он будет обозначаться одним символом.

Схема алгоритма решения системы однородных уравнений состояния показана на рис. 9.21. Исходными данными для алгоритма являются матрицы А, С и вектор начальных условий  Интегрирование уравнений состояния ведется с шагом .

Условный оператор конца счета должен формироваться исходя из решения конкретной задачи.

Итерационный алгоритм. Алгоритм (см. рис. 9.21) имеет один существенный недостаток — обычно шаг намного меньше, чем максимальное время исследования динамики Поэтому возникает необходимость многократного вычисления при достаточно больших  что требует значительного машинного времени. Имеется возможность сократить затраты машинного времени путем применения итерационного процесса вычисления когда для определения каждого последующего значения  используются результаты предыдущих вычислений. При этом достаточно вычислить матрицу перехода только для одного момента времени Сущность итерационного процесса заключается в следующем. Вектор переменных состояния рассчитанный для момента времени принимается за начальный вектор вместо и используется для расчета в момент времени т. е.

Повторяя этот процесс многократно, получаем последовательность значений вектора переменных состояния в моменты времени Эта последовательность в каждом цикле используется для вычисления

Схема итерационного алгоритма вычисления показана на рис. 9.22. Алгоритм вычисления находится вне цикла по и поэтому определение осуществляется один раз только для 

Такое построение алгоритма обусловливает накопление ошибок за счет неточного вычисления  Действительно, при использовании итерационной процедуры счета вектор переменных состояния при равен

Допустим, что матрица определена с матричной погрешностью  т. е.

где — точное значение матрицы Тогда

Ошибка  определяется выражением

Ошибка будет расти с увеличением

Итерационный алгоритм с компенсацией ошибок. Для компенсации накапливаемых ошибок в итерационном процессе можно применить метод уточнения вектора в момент С этой целью определяется не одна матрица перехода а две: При этом целесообразно применить соотношение

где — целое число.

Сначала итерационный процесс использует матрицу  для вычисления

При вектор переменных состояния вычисляется с использованием матрицы т. е.

При необходимости можно использовать не одну матрицу для уточнения вектора в моменты а несколько, значения которых рассчитаны для моментов времени Матрица компенсирует ошибки в моменты 

 матрица — в моменты . Точность вектора  определяется точностью матрицы и поэтому накопленные ошибки в соответствии с выражением (9.99) ликвидируются. В этот же момент времени

заменяется вектор начальных условий новое значение становится равным увеличивается на величину. Далее вычисление  снова происходит с матрицей до нового значения Когда текущее время t достигает значения  снова происходит компенсация ошибок.

Выбором соотношений между а также их количества можно обеспечить требуемую точность и малый шаг вычисления для достаточно длительных процессов при минимальной затрате машинного времени. На рис. 9.23 показана схема алгоритма, где применен однотактный способ компенсации накапливаемых ошибок. Сначала определяются матрицы и которые хранятся в памяти и используются для вычисления Два условных оператора  дают разветвление алгоритма, позволяющее использовать матрицу для уточнения вектора при последующим изменением момента

Построение алгоритмов решения неоднородных уравнений состояния

Алгоритмизации подлежит уравнение 

Первое его слагаемое отражает свободное движение системы за счет ненулевых начальных условий, второе и третье — определяют вынужденную составляющую решения.

Уравнение (9.103) удобно представить в виде следующих двух уравнений:

Уравнение (9.104) определяет вектор переменных состояния а (9.105) — выходной вектор по известным векторам и  Подставив  в (9.105), получим уравнение (9.103).

Итерационный алгоритм. Если принять то выражение (9.104) можно записать

где — шаг дискретности вычисления вектора

Для вычисления интеграла правой части уравнения (9.106) можно применять различные формулы приближенного вычислении определенных интегралов [4, 29]. Простой и удобной для применения является формула трапеции. Выбором шага дискретности можно обеспечить необходимую точность приближенного вычисления интеграла. Применение более сложных формул (Симпсона, Ньютона-Котеса, Гаусса, Маркова, Чебышева и др.) хотя и дает высокую точность интегрирования, однако затрудняет программирование и понимание построения алгоритма вычисления вектора состояния (9.106). Применение формулы трапеций позволяет выражения (9.106) и (9.105) привести к следующему виду:

Если принять за начальный вектор переменных состояния, то при равен

или после применения формулы трапеций

Сравнивая выражения (9.107) и (9.108) с (9.109) и (9.111), видим, что они идентичны по своему виду, за исключением аргументов вектора начальных условий и вектора выходных воздействий. Это дает возможность построить итерационный алгоритм вычисления переход к следующей итерации осуществляется изменением текущего времени на значение шага дискретности и вектора начальных условий. В каждой итерации используются только матрицы перехода при Поэтому достаточно вычислить эти матрицы раз.

Схема алгоритма показана на рис. 9.24.

Исходными данными для него являются матрицы А, В, С и D. вектор начальных условий начальное время  и шаг дискретности . Сначала рассчитывается матрица затем находятся которые используются для вычисления и в конечном итоге Логический оператор «анализ конца» проверяет условия выполнения конца расчета, заданного для конкретной задачи. При невыполнении этого условия присваивается новое значение вектору начальных условий и начинается следующая итерация счета. На входе можно иметь как отдельные составляющие выходного вектора, соответствующие слагаемым уравнения (9.106), так и суммарное значение вектора

Порядок исследования электрических цепей в среде Mathcad методом пространства состояния

(на примере схемы на рис. 925)

1. Составить согласно принципиальной схеме уравнение электрической цепи в дифференциальной форме.
2. Составить расширенную схему электрической цепи в переменных состояния.
3. Составить в стандартной форме уравнения состояния электрической цепи.
4. Экспериментально построить вектор выхода электрической цепи.

Исходные параметры системы уравнений состояния:

Параметры вектора входа:

Исходные параметры электрической цепи:

расширенной матрицы коэффициентов Ар — 
расширенной матрицы выхода Ср  —                 

вектора начальных условий —                              

Программа вычисления матрицы перехода

Примеры решения задач

Пример 9.9.1.

Построить методом пространства состояний переходные функции тока  напряжения на емкости и сопротивлении для электрической цепи, приведенной на рис. 9.25.

Дано:

Решение

1. Составляем уравнение электрической цепи

2. Составляем расширенную схему в переменных состояния для составления расширенных уравнений состояния. Уравнение (9.112) представим в виде

Расширенная схема в переменных состояния для решения уравнения (9.113), составленная по методу комбинирования производных, приведена на рис. 9.26.

По схеме в переменных состояния записываем уравнения состояния в виде

В качестве переменных вектора выхода принимаем напряжения и ток  Для этих переменных уравнения вектора выхода запишем в виде

В матричной форме уравнения состояния принимают вид

или — расширенная матрица коэффициентов электрической цепи.

В матричной форме уравнения вектора выхода принимают вид
      
или где — расширенная матрица выхода электрической цепи.

Вектор начальных условий по условию задачи принимает вид

Решение уравнений состояния и формирование вектора выхода проведем в среде Mathcad.

Исходные параметры системы

Параметры электрической цепи

Вводим символьное обозначение матриц:

Программа вычисления матрицы перехода

Матрица перехода при с принимает вид

Программа вычисления и построения вектора выхода

Переменные вектора выхода показаны па рис. 9.27.

Пример 9.9.2. 

Для электрический цепи (рис. 9.28) после замыкания ключа составить уравнение равновесия, расширенную схему в переменных состояния, уравнения состояния. Решить уравнения состояния в среде Mathcad методом пространства состояний. Вектор выхода сформировать из переменных 

Дано:

Решение

1. Составляем уравнение электрического равновесия для цепи в виде

или  

2.Составляем расширенную схему в переменных состояния. Уравнение равновесия преобразуем к виду

Схема в переменных состояния, составленная по методу комбинирования производных, представлена на рис. 9.29.

3. По виду схемы в переменных состояния составляем уравнения состояния. Схема в переменных состояния содержит два интегратора, следовательно, порядок уравнения На втором интеграторе реализована модель источника ЭДС 

Уравнения состояния в матричной форме принимают вид 

Расширенная матрица коэффициентов представляется в виде

4. Уравнения вектора выхода системы в матричной форме имеют вид

Расширенная матрица выхода и вектор начальных условий по условию задачи представлены соответственно в виде:

В среде Mathcad решаем систему уравнений состояния и из переменных состояния формируем вектор выхода электрической цепи.

Выходные переменные приведены на рис. 9.30.

Пример 9.9.3.

Для электрической цепи, приведенной на рис.9.31, после замыкания ключа составить уравнение равновесия, расширенную схему в переменных состояния. Решить уравнения состояния в среде Mathcad. Вектор выхода сформировать из переменных состояния:

Дано:

Решение

1.Составляем по второму закону Кирхгофа уравнение электрического равновесия для цепи после ее переключения в виде

Так как  то после подстановки уравнение примет вид

2. Составляем схему электрической цепи в переменных состояния. Так как после переключения внешний источник ЭДС отсутствует, т.е. то расширения схемы не требуется.

Разрешим уравнение равновесия относительно старшей производной, тогда

а схема в переменных состояния для решения этого уравнения имеет вид (рис. 9.32).

3. По схеме в переменных состояния (рис. 9.32) составляем уравнения состояния:

или в матричной форме

Матриц коэффициентов определяется выражением 

4. Формируем из переменных состояния вектор выхода. Уравнения выхода имеют вид

В матричной форме

Матрица выхода и вектор начальных условий принимают вид:

5.Решение уравнений состояния можно выполнить любым численным методом. На рис. 9.33 приведены переменные вектора выхода, вычисленные методом пространства состояний в среде Mathcad по программе, приведенной в задаче 9.9.1.

Пример 9.9.4.

Задана электрическая цепь (рис.9.34). Требуется составить уравнение равновесия цепи по второму закону Кирхгофа после замыкания ключа и расширенную схему в переменных состояния. Получить и решить уравнения состояния любым численным методом с применением ПЭВМ, сформировать вектор выхода из переменных состояния:

 и построить их графики.

 Дано: 

Решение

1. Записываем уравнение электрического баланса цепи в виде 

2. Составляем расширенную схему в переменных состояния при Уравнение равновесия запишем в виде

С учетом входного воздействия схема в переменных состояния электрической цепи может быть представлена в виде, приведенном на рис. 9.35.

3. В соответствии со схемой в переменных состояния запишем уравнении состояния электрической цепи

В матричной форме уравнения состояния имеют вид

Расширенная матрица коэффициентов представляется в виде

Вектор начального состояния из условия задачи запишем в виде
 

4. Определяем уравнения выхода из заданных переменных состояния электрической цепи.

В матричной форме уравнения вектора выхода имеют вид
  
Матрица коэффициентов выхода

5. Решение уравнений состояния выполняем в среде Mathcad с помощью программы из задачи 9.9.1.

Переменные вектора выхода приведены на рис. 9.36.

Из переменных вектора выхода наглядно видно, что реакция последовательного колебательного контура на единичное воздействие — переходная функция — имеет слабозатухающий колебательный характер.

Пример 9.9.5.

По электрической схеме, представленной на рис.9.37, после замыкания ключа составить уравнение электрического равновесия по второму закону Кирхгофа, расширенную схему в переменных состояния, уравнения состояния и выхода из переменных состояния

Дано:

Решить уравнения состояния одним из численных методов в среде Mathcad и построить графики компонентов вектора выхода.

Решение

1. Записываем уравнение электрического равновесия для электрической цепи

2.Составляем расширению схему в переменных состояния при  Схема приведена на рис. 9.38.

Схема модели гармонического воздействия выполнена на интеграторах 3 и 4, а модель цепи реализована на интеграторах 1 и 2.

3. Составляем уравнения состояния по схеме в переменных состояния

Запишем уравнения и матричной форме в виде 

Расширенная матрица коэффициентов системы имеет вид 

Вектор начального состояния системы для

4. Составляем уравнения вектора выхода из переменных состояния

В матричной форме уравнения выхода имеют вид

Расширенная матрица выхода запишется и виде

5. Решение уравнений состояния выполним в среде Mathcad по программе, представленной в задаче 9.9.1. Переменные вектора выхода представлены на рис. 9.39.

(Андрей Андреевич Марков (1856-1922) – русский математик, академик)

Определение. Процесс, протекающий в физической системе, называется Марковским, если в любой момент времени вероятность любого состояния системы в будущем зависит только от состояния системы в текущий момент и не зависит от того, каким образом система пришла в это состояние.

Определение. Цепью Маркова называется последовательность испытаний, в каждом из которых появляется только одно из K несовместных событий Ai из полной группы. При этом условная вероятность Pij(S) того, что в S –ом испытании наступит событие Aj при условии, что в (S – 1) – ом испытании наступило событие Ai, не зависит от результатов предшествующих испытаний.

Независимые испытания являются частным случаем цепи Маркова. События называются Состояниями системы, а испытания – Изменениями состояний системы.

По характеру изменений состояний цепи Маркова можно разделить на две группы.

Определение. Цепью Маркова с дискретным временем Называется цепь, изменение состояний которой происходит в определенные фиксированные моменты времени. Цепью Маркова с непрерывным временем Называется цепь, изменение состояний которой возможно в любые случайные моменты времени.

Определение. Однородной Называется цепь Маркова, если условная вероятность Pij перехода системы из состояния I В состояние J не зависит от номера испытания. Вероятность Pij называется Переходной вероятностью.

Допустим, число состояний конечно и равно K.

Тогда матрица, составленная из условных вероятностей перехода будет иметь вид:

Эта матрица называется Матрицей перехода системы.

Т. к. в каждой строке содержаться вероятности событий, которые образуют полную группу, то, очевидно, что сумма элементов каждой строки матрицы равна единице.

На основе матрицы перехода системы можно построить так называемый Граф состояний системы, его еще называют Размеченный граф состояний. Это удобно для наглядного представления цепи. Порядок построения граф рассмотрим на примере.

Пример. По заданной матрице перехода построить граф состояний.

Т. к. матрица четвертого порядка, то, соответственно, система имеет 4 возможных состояния.

S1

0,2 0,7

S2 0,4 S4

0,6 0,5

0,1 0,5

S3

На графе не отмечаются вероятности перехода системы из одного состояния в то же самое. При рассмотрении конкретных систем удобно сначала построить граф состояний, затем определить вероятность переходов системы из одного состояния в то же самое (исходя из требования равенства единице суммы элементов строк матрицы), а потом составить матрицу переходов системы.

Пусть Pij(N) – вероятность того, что в результате N испытаний система перейдет из состояния I в состояние J, R – некоторое промежуточное состояние между состояниями I И J. Вероятности перехода из одного состояния в другое Pij(1) = Pij.

Тогда вероятность Pij(N) может быть найдена по формуле, называемой Равенством Маркова:

Здесь Т – число шагов (испытаний), за которое система перешла из состояния I В состояние R.

В принципе, равенство Маркова есть ни что иное как несколько видоизменная формула полной вероятности.

Зная переходные вероятности (т. е. зная матрицу перехода Р1), можно найти вероятности перехода из состояния в состояние за два шага Pij(2), т. е. матрицу Р2, зная ее – найти матрицу Р3, и т. д.

Непосредственное применений полученной выше формулы не очень удобно, поэтому, можно воспользоваться приемами матричного исчисления (ведь эта формула по сути – не что иное как формула перемножения двух матриц).

Тогда в общем виде можно записать:

Вообще то этот факт обычно формулируется в виде теоремы, однако, ее доказательство достаточно простое, поэтому приводить его не буду.

Пример. Задана матрица переходов Р1. Найти матрицу Р3.

Определение. Матрицы, суммы элементов всех строк которых равны единице, называются Стохастическими. Если при некотором П все элементы матрицы Рп не равны нулю, то такая матрица переходов называется Регулярной.

Другими словами, регулярные матрицы переходов задают цепь Маркова, в которой каждое состояние может быть достигнуто через П шагов из любого состояния. Такие цепи Маркова также называются Регулярными.

Теорема. (теорема о предельных вероятностях) Пусть дана регулярная цепь Маркова с п состояниями и Р – ее матрица вероятностей перехода. Тогда существует предел и матрица Р(¥) имеет вид:

Т. е. матрица состоит из одинаковых строк.

Теперь о величинах Ui. Числа U1, U2, …, Un называются Предельными вероятностями. Эти вероятности не зависят от исходного состояния системы и являются компонентами собственного вектора матрицы РТ (транспонированной к матрице Р).

Этот вектор полностью определяется из условий:

Пример. Найдем предельные вероятности для рассмотренного Выше примера.

C учетом того, что U1 + U2 = 1, получаем:

Получаем:

Соседние файлы в предмете [НЕСОРТИРОВАННОЕ]

• #
• #
• #
• #
• #
• #
• #
• #
• #
• #
• #