Идентификация и моделирование систем управления (стр. 2 из 3)
Напряжение на участке исключенного источника тока uL =Е, т. к. тока в разомкнутой цепи нет, то следовательно нет падения напряжения на активных сопротивлениях.
После получения всех коэффициентов матриц A и B можно записать систему (2. 2) для полученных коэффициентов: a11 = (- 1 / C )(1/ (R2 + R3)); a12 = (1 / C )(R2 / (R2 + R3)); a21 = -R2 /L(R2 + R3); a22 = -(1/L)( R1 + (R2 R3 / R2 + R3)); b11 =0; b21 = 1/L.
Подставляя в полученную систему численные значения параметров компонент, согласно исходной схеме, получим.
В координатной форме полученная система имеет вид
Возвращаясь к первоначальным переменным x1 = uc; x2 = iL , можно записать в общем виде для заданной электрической цепи следующую систему уравнений в форме Коши, которую необходимо решить и выполнить анализ динамического процесса с помощью средств система автоматизации математических расчетов MATLAB и пакета динамических систем Simulink, входящего в состав расширенных версий MATLAB, а также вручную и сравнить полученные результаты, которые должны совпасть.
Система (2.3)
2.2 Анализ динамических процессов в системе на основе использования построенной аналитической модели
После получения динамической модели изучаемой системы в виде системы дифференциальных уравнений, записанных в форме пространства состояний, необходимо выполнить анализ динамических процессов протекающих в системе. Для выполнения этой задачи следует найти решение системы уравнений, т.е. найти аналитическое выражение – функцию, отражающую закон, согласно которому изменяются переменные состояния во времени. Получив закон, можно определить характер динамических процессов, протекающих в системе.
 |
Для нахождения решения системы (2.3):
используем матрично-векторное соотношение
x(t)= А-1 (exp(Аt) - 1 ) f0 , (2.4)
имеющее место при нулевых начальных условиях и внешнем воздействии f0 в виде вектора с постоянными компонентами.
Сначала находим А-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
Полученный полином будет определяться выражением
λ2 + 169.79λ + 315.
Соответственно, его корни:
>> A=[-1.79 7.14;-2 -168];
>> [D]=eig(A)
D = -1.8760
-167.9140
Следовательно: λ1=--1.8760; λ2 = -167.914
Результат вычислений позволяет сделать вывод о апериодическом характере переходного процесса и устойчивости системы, необходимым и достаточным условием которой является отрицательность корней характеристического уравнения.
Для вычисления матричной экспоненты используется формула Сильвестра, согласно которой
Для системы (2.4) можно записать экспоненту:
Подставляя все полученные данные в соотношение (2.4) приходим к следующему выражению:
После выполнения действий над матрицами получим следующее решение системы (1):
x1(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
x2(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
Для построения графиков функций x1(t) и x2(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 (изменения во времени переменных состояния системы x1(t) и x2(t) носят апериодический характер). На рис. 3 приведена распечатка графика функции
x1(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
при начальном условии x1 (0) = 0 (график начинается с ординаты x1 = 0).
На рис. 4 приведена распечатка графика функции
x2(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
при начальном условии x2 (0) = 0 (график начинается с ординаты x2 = 0).
Полученные графики демонстрируют апериодический переходный процесс, возникающий в электрической цепи при подключении источника постоянной э.д.с. При этом напряжение на конденсаторе x1 = uc и ток через катушку индуктивности x2 = iL изменяются согласно найденных соотношений:
x1(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
и
x2(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,
соответственно.
Как видно из графиков временных зависимостей, процесс асимптотически приближается к установившемуся состоянию с принужденной составляющей x1(t)= x1 =-4.7 и x2(t)= x2 = -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 будет иметь следующий вид. Начальные условия: x1 (0) = 0; x2 (0) = 0.
>> [t, x]=ode45('difur',[0 10],[0;0])
Решение исходной системы (2.3) в виде числовых массивов сохраняются в папке work в виде бинарного файла с расширением .mat. Для этого служит команда Save Workspace As…. в меню File. При необходимости сохраненный файл можно загрузить в рабочую область Workspace (командой load).