СОДЕРЖАНИЕ

1. Задание           3

2. Построение  аналитической модели и ее  анализ.

     2.1 Построение аналитической модели    4

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

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

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

     2.5. Моделирование с использованием имитационного пакета  моделирования динамических систем Simulink     25 
 
 
 
 
 
 
 
 
 
 
 
 

1. Задание

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

     Схема электрической цепи и параметры составляющих ее компонент:

                                                  Рис. 1. Электрическая RLC - цепь 
 
 
 
 
 
 
 
 
 

     2. Построение аналитической  модели и ее анализ.

     2.1 Построение аналитической модели 

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

            Известно несколько способов  составления уравнений состояния.  Рассмотрим наиболее целесообразный  способ, основанный на сведении  послекоммутационной схемы к  резистивной с источниками э.д.с. и тока. С этой целью индуктивные элементы в послекоммутационной схеме заменяют на источники тока, которые возбуждают ток в том же направлении, что и в исходной схеме, а конденсатор на источник э.д.с. с э.д.с. направленной встречно току в ветви с конденсатором, т. е. встречно u_c . В результате схема окажется без реактивных элементов (резистивной), но с дополнительными источниками тока и э.д.с.

           В полученной резистивной схеме  один из узлов заземляют и  составляют уравнения по методу  узловых потенциалов. По найденным потенциалам узлов рассчитывают напряжения на источниках тока, эквивалентных индуктивным элементам и токи через источники э.д.с., эквивалентные емкостным элементам. Далее разрешают уравнения цепи относительно производных di_L /dt и du_C /dt и получают запись системы дифференциальных уравнений в нормальной форме (форма Коши).

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

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

           Рассмотрим модель электрической цепи (Рис. 1).

           Для построения модели выбираем  вектор переменных состояния 

     т. е. x₁ = u_c ; x₂ = i_L. Так как u_L = L(di_L /dt), то dx₂/dt = di_L /dt = u_L /L, а u_c = 1/C i_c dt, то dx₁/dt = du_c /dt = i_c / C.

     Уравнения состояния в матричной форме  в общем виде для 

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

      ( 2.1)

           Коэффициенты матриц будем определять методом наложения при рассмотрении эквивалентной резистивной схемы (Рис. 2).

                              Рис. 2. Эквивалентная резистивная схема. 

     Запишем систему (2.1) в координатной форме, из которой определим коэффициенты матриц A и B.

                         Система (2.2) 

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

     моделирования такого режима необходимо в эквивалентной  резистивной схеме закоротить источник э.д.с. u₁ = E = 0. А для определения коэффициента  a₁₁ матрицы A исключить источник тока x₂ = i_L = 0. Соответственно из первого уравнения системы (2. 2) получим:

     При условии:  E = 0; i_L = 0.

     Для измененной схемы определяем i_c = - u_c / (R₂ + R₃) и подставляем в выражение для a₁₁ .В результате получим:

     Для определения коэффициента  a₁₂ матрицы A восстанавливаем источник тока, но замыкаем источник э. д. с.: x₁ = u_c = 0. Соответственно из первого уравнения системы    (2. 2) получим:

     При условии:  E = 0; u_c  = 0.

            Для измененной схемы определяем i_c =(R₂ /(R₂ + R₃)) i_L

       и подставляем в выражение для a₁₂ .В результате получим:

     Для определения коэффициента a₂₁ матрицы A в эквивалентной резистивной схеме закоротить источник э.д.с. u₁ = E = 0 и исключаем источник тока x₂ = i_L = 0.

          Соответственно из второго уравнения системы получим:

     При условии:  E = 0; i_L = 0.

           Для измененной схемы определяем  u_L = i_c R₂ , для этого находим i_c = - u_c / (R₂ + R₃) и подставляем в u_L = i_c R₂ . Получим: u_L =( - u_c / (R₂ + R₃))R₂ . Тогда

     Для определения коэффициента a₂₂ матрицы A в эквивалентной резистивной схеме закоротить источник э.д.с. u₁ = E = 0, восстанавливаем источник тока, но замыкаем источник э. д. с.: x₁ = u_c = 0. Соответственно из второго уравнения системы получим:

     Для измененной схемы определяем: u_L = - (R₁ + (R₂ R₃ / R₂ + R₃)) i_L . Тогда a₂₂ = -(1/L)( R₁ + (R₂ R₃ / R₂ + R₃)).

           Определение коэффициентов матрицы  B (коэффициенты матрицы B определяют вклад входных величин в баланс токов и напряжений) предполагает исключение источника тока x₂ = i_L = 0, замыкание источника э. д. с. x₁ = u_c = 0 и сохранение источника u₁ = E. Тогда для определения коэффициента b₁₁ матрицы B в первом уравнении системы  полагаем x₁ = u_c = 0, x₂ = i_L = 0. Получим:

но i_c =0 при любом Е, т. к. ветвь с источником тока разомкнута, то b₁₁ =0.

     Для определения коэффициента b₂₁ матрицы B во втором  уравнении системы (4. 2) полагаем x₁ = u_c = 0, x₂ = i_L = 0, что предполагает исключение источника тока x₂ = i_L = 0, замыкание источника э. д. с. x₁ = u_c = 0 и сохранение источника u₁ = E.  Получим:

     Напряжение  на участке исключенного источника  тока  u_L =Е, т. к. тока в разомкнутой цепи нет, то следовательно нет падения напряжения на активных сопротивлениях.

     После получения всех коэффициентов матриц  A и B можно записать систему (2. 2) для полученных коэффициентов: a₁₁ = (- 1 / C )(1/ (R₂ + R₃)); a₁₂ = (1 / C )(R₂ / (R₂ + R₃)); a₂₁ = -R₂ /L(R₂ + R₃); a₂₂ = -(1/L)( R₁ + (R₂ R₃ / R₂ + R₃)); b₁₁ =0; b₂₁ = 1/L.

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

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

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

           Система (2.3) 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

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

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

      Для нахождения решения системы (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):