Идентификация и моделирование систем управления

      Для нахождения решения системы (2.3):

     

       

     используем  матрично-векторное соотношение

                              x(t)=  А^-1 (exp(Аt) - 1 ) f₀ ,   (2.4)

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

       Сначала находим А^-1 по известной  А  (в МatLab).

>> A=[-1.79 7.14;-2 -168];

>> inv(A)

ans = -0.5333   -0.0227

            0.0063   -0.0057

      Далее находим матричную экспоненту  exp(Аt), составив предварительно характеристическое уравнение системы   (2.3) и найдя его корни

>> poly(A)

ans = 1.0000  169.7900  315.0000

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

                                           λ² + 169.79λ + 315.

     Соответственно, его корни:

     >> A=[-1.79 7.14;-2 -168];

     >> [D]=eig(A)

     D =  -1.8760

             -167.9140

     Следовательно: λ₁=--1.8760; λ₂ = -167.914

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

           Для вычисления матричной экспоненты используется формула Сильвестра, согласно которой

       
 
 
 

     Для системы (2.4) можно записать экспоненту:

     

     Подставляя  все полученные данные в соотношение (2.4) приходим к следующему выражению:

     

     После выполнения действий над матрицами получим следующее решение системы (1):

     x₁(t)=-16621/400*(-3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200+1/200*275686441^(1/2))*t)-1/400*(-3260759/82705932300*275686441^(1/2)+179/300)*275686441^(1/2)*exp((-16979/200+1/200*275686441^(1/2))*t)-16621/400*(3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200-1/200*275686441^(1/2))*t)+1/400*(3260759/82705932300*275686441^(1/2)+179/300)*275686441^(1/2)*exp((-16979/200-1/200*275686441^(1/2))*t)-119/25

     x₂(t)=(-3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200+1/200*275686441^(1/2))*t)+(3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200-1/200*275686441^(1/2))*t)-179/150

     Для построения графиков функций x₁(t) и x₂(t) выполним команду plot в режиме командной строки. Программы приведены ниже.

     >> t=0:0.1:12.54;

     >>x=-16621/400*(-3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200+1/200*275686441^(1/2))*t)-1/400*(-3260759/82705932300*275686441^(1/2)+179/300)*275686441^(1/2)*exp((-16979/200+1/200*275686441^(1/2))*t)-16621/400*(3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200-1/200*275686441^(1/2))*t)+1/400*(3260759/82705932300*275686441^(1/2)+179/300)*275686441^(1/2)*exp((-16979/200-1/200*275686441^(1/2))*t)-119/25;

     >> plot(t,x)

     

     Рис 3

     >> t=0:0.02:5;

     >>x=-(-3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200+1/200*275686441^(1/2))*t)+(3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200-1/200*275686441^(1/2))*t)-179/150;

     >> plot(t,x)

     

     Рис 4

           Результат моделирования, отражающий динамический процесс представлен на рис. 3 и рис. 4 (изменения во времени переменных состояния системы x₁(t) и x₂(t) носят апериодический характер). На рис. 3  приведена распечатка графика функции

     x₁(t)=-16621/400*(-3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200+1/200*275686441^(1/2))*t)-1/400*(-3260759/82705932300*275686441^(1/2)+179/300)*275686441^(1/2)*exp((-16979/200+1/200*275686441^(1/2))*t)-16621/400*(3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200-1/200*275686441^(1/2))*t)+1/400*(3260759/82705932300*275686441^(1/2)+179/300)*275686441^(1/2)*exp((-16979/200-1/200*275686441^(1/2))*t)-119/25

 при начальном  условии  x₁ (0) = 0 (график начинается с ординаты  x₁ = 0).

На рис. 4 приведена распечатка графика функции

       x₂(t)=(-3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200+1/200*275686441^(1/2))*t)+(3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200-1/200*275686441^(1/2))*t)-179/150

при начальном  условии x₂ (0) = 0 (график начинается с ординаты  x₂ = 0).

     Полученные  графики демонстрируют апериодический переходный процесс, возникающий в  электрической цепи при подключении  источника постоянной э.д.с. При этом напряжение на конденсаторе x₁ = u_c и ток через катушку индуктивности x₂ = i_L  изменяются согласно найденных соотношений:

     x₁(t)=-16621/400*(-3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200+1/200*275686441^(1/2))*t)-1/400*(-3260759/82705932300*275686441^(1/2)+179/300)*275686441^(1/2)*exp((-16979/200+1/200*275686441^(1/2))*t)-16621/400*(3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200-1/200*275686441^(1/2))*t)+1/400*(3260759/82705932300*275686441^(1/2)+179/300)*275686441^(1/2)*exp((-16979/200-1/200*275686441^(1/2))*t)-119/25

     и

     x₂(t)=(-3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200+1/200*275686441^(1/2))*t)+(3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200-1/200*275686441^(1/2))*t)-179/150,

       соответственно.

           Как видно из графиков временных  зависимостей, процесс асимптотически приближается к установившемуся состоянию с принужденной составляющей x₁(t)= x₁ =-4.7 и x₂(t)= x₂ = -1.2.

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

    

     2.3. Моделирование с использованием солверов

   Следующий этап моделирования состоит в  нахождении решений исходной системы  (2.3), используя встроенные средства для решения дифференциальных уравнений и их систем. Система MATLAB выполняет численное решение обыкновенных дифференциальных уравнений произвольного порядка и систем с начальными условиями. В MATLAB имеется целый ряд встроенных функций, предназначенных для решения задачи Коши для обыкновенных дифференциальных уравнений. Библиотека включает несколько функций, реализующих различные методы решения задачи Коши ode (ordinary differential equations).

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

            Для  системы (2. 3) напишем текст программы, для этого в меню File окна системы  MATLAB выполним команду New m-file и в открывшемся окне редактора/отладчика m-файлов наберем текст файл-функции, который будет таким.

                         function F=difur(t,x)

                         F=[-1.79*x(1)+7.14*x(2);-2*x(1)-168*x(2)-210];

           Сохраним его в файле difur.m в текущем каталоге.

           Для вычисления решения системы  на интервале [0, 10] используем  командную строку. С учетом начальных условий, обращение к функции ode 45 будет иметь следующий вид. Начальные условия: x₁ (0) = 0;  x₂ (0) = 0.

                    >> [t, x]=ode45('difur',[0 10],[0;0])

   Решение исходной системы (2.3) в виде числовых массивов сохраняются в папке work в виде бинарного файла с расширением .mat. Для этого служит команда Save Workspace As…. в меню File. При необходимости  сохраненный файл можно загрузить в рабочую область Workspace (командой load).

   После выполнения сеанса работы (сессии) в  окне команд Command Window появятся результаты решения системы (1) в виде массивов аргумента t и искомых функций x, которые можно сохранить в папке work под именем rech1. mat.

           Для отображения графика исходной  системы дифференциальных уравнений (2.3) необходимо в режиме командной строки выполнить команду plot:

                             >> plot(t,x(:,1),'r',t,x(:,2),'k--')

     

     Рис 5

           На рис. 5 приведена экранная форма графика решения заданной системы уравнений (2. 3) x₁(t), x₂(t).

           Графики интерпретируют изменение  во времени величины функции: 

     x₁(t)=-16621/400*(-3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200+1/200*275686441^(1/2))*t)-1/400*(-3260759/82705932300*275686441^(1/2)+179/300)*275686441^(1/2)*exp((-16979/200+1/200*275686441^(1/2))*t)-16621/400*(3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200-1/200*275686441^(1/2))*t)+1/400*(3260759/82705932300*275686441^(1/2)+179/300)*275686441^(1/2)*exp((-16979/200-1/200*275686441^(1/2))*t)-119/25 и

     x₂(t)=(-3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200+1/200*275686441^(1/2))*t)+(3260759/82705932300*275686441^(1/2)+179/300)*exp((-16979/200-1/200*275686441^(1/2))*t)-179/150, соответственно.

           Как видно из графиков временных  зависимостей, процесс асимптотически  приближается к установившемуся  состоянию с принужденной составляющей  x₁(t)= x₁ =-4.7 и x₂(t)= x₂ = -1.2.

     Сравнение полученных графиков с ранее построенными x₁(t) x₂(t) (Рис. 5 ) при моделировании в режиме командной строки свидетельствует о идентичности результатов, полученных при использовании различных способов взаимодействия с программой при реализации ее широких возможностей.

2.4. Моделирование с использованием пакета расширения Symbolic Math Tolbox

           Возможности пакета расширения  Symbolic Math Tolbox позволяют в рамках системы MATLAB осуществлять аналитические вычисления и аналитические преобразования выражений. Помимо выполнения аналитических преобразований, пакет Symbolic Math Tolbox позволяет выполнять арифметические вычисления с контролируемой точностью, которую можно задать заранее. Этот пакет осуществляет алгебраические операции над объектами нового типа – sym – объектами. Такие объекты получаются после вызова одноименной функции конструктора таких объектов. Над объектами типа sym производятся манипуляции в соответствии с правилами алгебры и математического анализа.

   Для решения дифференциальных уравнений  в форме Коши MATLAB имеет следующую функцию:

   dsolve(‘eqn1’,’eqn2’,….) – возвращает аналитическое решение системы дифференциальных уравнений с начальными условиями. Сначала задаются уравнения, затем начальные условия (равенствами eqn).

   По  умолчанию в качестве независимой  переменной задается переменная ‘t’. Можно использовать и другую переменную, добавив ее в конце списка параметров функции dsolve. Символ D обозначает производную по независимой переменной, D2 означает вторую производную и т. д.

   Начальные условия задаются в виде равенств ‘y(a)=b’, ‘Dy(a)=b’ , где y – независимая переменная, a и b – константы. Если число начальных условий меньше, чем число дифференциальных уравнений, то в решении будут присутствовать произвольные постоянные С1, С2,…и т. д.

   Выполним  функцию dsolve для системы (1):