Использование ============= .. _installation: Установка --------- В этой версии нужно иметь следующие файлы (в каталоге, куда настроены пути импорта): * libKernel.py * libComm.py * libGraph2D.py Зависимости ~~~~~~~~~~~ * matplotlib * numpy * pyquaternion Для запуска примеров: * scipy * control * cobs * pyserial Установите их своим пакетным менеджером (mamba, conda или pip): .. code-block:: console (.venv) $ mamba install -c conda-forge matplotlib numpy pyquaternion scipy control slycot pyserial git pip (.venv) $ pip install git+git://github.com/cmcqueen/cobs-python.git .. _usage: Примеры ------- Рассмотрим несколько вариантов моделей и сравним их с аналогичными моделями, выполненными в SimInTech и с теоретически рассчитанными переходными процессами. Передаточная функция ~~~~~~~~~~~~~~~~~~~~ Смоделируем переходной процесс в передаточной функции (ПФ) при нулевых начальных условиях от входного воздействия типа единичный скачок. Для понимания кода нужно иметь опыт работы с библиотекой `pycontrol `_. Используем ПФ вида .. math:: W(s) = \dfrac{s+1}{0.01s^2 + 0.01s + 1}. Обозначим входное воздействие переменной u, а выход - переменной y. База данных будет выглядеть следующим образом: :: class DataBase(): # --- общие переменные --- t = 0. # время dt = 0.25 # шаг задач tmax = 4. # время моделирования cmd = 0 # команда всем задачам y = 0.0 # переходной процесс u = 1.0 # входное воздействие Словарь Tasks с описанием задач и переменных для обмена выглядит следующим образом (здесь только наша единственная задача TF - имя совпадает с именем класса): :: Tasks = { 'TF': {'Keys':'t,cmd,dt,y,u'}, } Чтобы модель могла обмениваться данными с диспетчером, необходимо её запустить - следующими двумя строчками после определения класса модели (здесь бд и список задач определены в том же файле, что и задача): :: model = TF(TaskList=Tasks, DB=DataBase(), isSheduler=True) model.Manager.Loop() Рассчитаем переходной процесс двумя способами: 1. Разобьём наш интервал времени на шаги равные dt, и будем считать решение только в этих точках. 2. Каждый из шагов разобьём на микрошаги, чтобы получить более детальный переходной процесс. .. Это достаточно учебный пример, для демонстрации способов получения решения дифференциальных уравнений Способ с большими шагами ************************ Код задачи для первого способа будет выглядеть следующим образом (пример ex1): :: class TF(libKernel.TaskTemplate): """ объект в форме передаточной функции """ def Setup(self): self.Num = [80] self.Den = [1, 0, 0] self.System = ct.tf2ss(self.Num, self.Den) self.plot1 = Plot2D(DB=self) self.plot1.Setup('t', [['y']) def Initialize(self): # задаём начальное состояние (нулевые НУ, но можно задавать из БД) self.State = np.zeros((len(self.Den)-1, 1)) self.u = 1.0 # управление =const self.plot1.Initialize() # инициализируем график # рассчитаем начальное значение на выходе res = ct.initial_response(self.System, T=[0], X0=self.State, T_num=2) self.y = res.outputs[0] self.plot1.Run() # добавляем первую точку на график def Run(self): # будем считать по точкам, на интервал dt _, self.y, self.State = ct.forced_response(self.System, T=[0, self.dt], U=self.u*np.ones(2), return_x=True, squeeze=True, X0=self.State[:,-1]) self.y = self.y[-1] self.plot1.Run() def Finalize(self): self.plot1.Finalize() Для расчёта переходного процесса используется функция `forced_response `_, вызов которой находится самом начале метода Run. Она возвращает переходной процесс (y) и cостояние системы (State) в последний момент интервала времени. Моменты времени, в которых будет рассчитано решение, задаётся параметром Т - вектором из двух точек [0, dt], времен начала и конца расчётного интервала. Первым параметром для forced_response указывается система System, в форме пространства состояний (forced_response эффективнее для такой формы). Входное воздействие задано переменной u = 1.0, в методе Initialize. В Initialize удобно вынести параметры задачи, не меняющиеся при расчёте. Там же задано начальное состояние - вектор State, размерность которого равна порядку уравнения (или количеству интеграторов системы, что то же самое). В методе Setup задаётся система System, с числителем Num и знаменателем Den (используем функцию `tf2ss `_, чтобы получить ПФ в пространстве состояний). В этой задаче для отображения переходного процесса использован двумерный "бегущий" график (отображает переходной процесс по мере его расчёта). Он инициализируется в методе Setup, где создаётся объект plot1, экземпляр класса Plot2D. При создании указывается база данных DB, из которой берутся переменные для отображения на графике (t и y). В качестве базы указывается собственно класс нашей задачи - self (именно в нём менеджер данных заводит атрибуты - t и y). В методе Initialize и методе Run вызывается метод Plot2D.Run(). Он автоматически берет значение указанной (при настройке) переменной и добавляет её на график - точку НУ (начальный элемент вектора результата) или на конце интервала dt. В методе Finalize график финализируется - при окончании расчёта задачи график дорисовывает все точки, которые он успел получить. В итоге откроется такой график: .. image:: img/tTFR.png :alt: Переходной процесс для интервалов dt=0.1c. Микрошаги ********* Это способ с промежуточными шагами (пример ex2.1). Код задачи отличается только параметрами расчёта и отрисовки: :: class TF(libKernel.TaskTemplate): """ объект в форме передаточной функции, микрошаги """ def Setup(self): self.Np = 50 # количество микрошагов в интервале dt self.Num = [1, 1] self.Den = [0.01, 0.01, 1] self.System = ct.tf2ss(self.Num, self.Den) # вектор времени self.tSol = np.linspace(0, self.dt, self.Np) self.u = np.ones(self.tSol.shape) # управление =const self.plot1 = Plot2D(DB=self) self.plot1.Setup('time', [['y']], chunkLen=self.Np) def Initialize(self): # задаём ПФ и начальное состояние (нулевые НУ) self.t0 = 0. self.State = np.zeros((len(self.Den)-1, 1)) self.plot1.Initialize() def Run(self): _, self.y, self.State = ct.forced_response(self.System, T=self.tSol, U=self.u, return_x=True, squeeze=True, X0=self.State[:,-1]) self.time = self.tSol + self.t0 self.t0 += self.dt self.plot1.Run() self.y = self.y[-1] # в БД пишем решение в конце интервала def Finalize(self): self.plot1.Finalize() В методе Setup задаём количество микрошагов (Np=50) в интервале dt, и вектор времени tSol. И поскольку мы знаем, что входное воздействие (скачок) не будет меняться при моделировании, зададим его сразу - вектором с единичками, длина которого совпадает с tSol. Обратите внимание, что настройка графика усложняется. Задаётся ещё один параметр - длина чанка chunkLen, то есть количество точек, которые будут поступать на график каждый расчётный шаг. Остальной код примерно такой же. Единственное, что в функцию расчёта forced_response передаётся вектор времени tSol. В последней строчке метода Run мы записываем в БД последнюю точку решения - в данном случае база настроена на обмен единственным значением y на каждый шаг dt. Полное решение отображается только на графике. Его можно сохранить для дальнейшего анализа, если нажать на кнопку Save. В итоге откроется график переходного процессп для тех же интервалов по 0.5c, но с микрошагами: .. image:: img/tTF.png :alt: Переходной процесс для интервалов dt=0.5c с микрошагами. Совместив оба графика, видим, что точки решения совпадают в моменты времени, кратные dt: .. image:: img/tTFres.png :alt: Расчёты совпадают. Значение выходной переменной y и при большом шаге, и при разбивке на микрошаги получилось одинаковым - потому что библиотека control для систем с непрерывным временем использует аналитический расчёт ПП `с использованием матричной экспоненты `_. Дифференциальное уравнение ~~~~~~~~~~~~~~~~~~~~~~~~~~ В случае, когда требуется решить дифференциальное уравнение общего вида, удобно применить функцию `odeint `_ библиотеки scipy. Наша передаточная функция в форме Коши (или в форме пространства состояний) будет иметь вид .. math:: \begin{eqnarray} \dot{X} &=& A \cdot X + B \cdot U \\ Y & = & C \cdot X + D \cdot U, \end{eqnarray} где .. math:: \begin{alignat}{2} A & = \left[\begin{array}{cc} -1 & 10\\ -10& 0 \end{array} \right] & \qquad B & = \left[\begin{array}{c} 1\\ 0 \end{array} \right] \\ C & = \left[\begin{array}{cc} 100 & -10 \end{array} \right] & D & = \left[ 0 \right]. \end{alignat} Зададим правую часть ДУ: :: A = np.array([[-1., 10.], [-10.0, 0]]) B = np.array([[1.], [0.]]).T C = np.array([[100., -10.]]) D = np.array([[0.]]) def de(y, t): dydt = np.matmul(A, y) + B return np.squeeze(dydt) .. Обратим внимание, что что матрица B транспонирована, а в функции правых частей de применяется Код задачи выглядит следующим образом (пример ex2.2): :: class ODE(libKernel.TaskTemplate): """ используем для решения scipy.integrate.odeint """ def Setup(self): self.Np = 50 self.tSol = np.linspace(0, self.dt, self.Np) # вектор времени self.plot1 = Plot2D(DB=self) self.plot1.Setup('time', [['y']], chunkLen=self.Np) def Initialize(self): # задаём начальное состояние (нулевые НУ) self.t0 = 0. self.State = [0, 0] self.plot1.Initialize() def Run(self): sol = odeint(de, self.State, self.tSol) Y = np.matmul(C, sol.T) + D self.y = np.squeeze(Y) self.time = self.tSol + self.t0 self.t0 += self.dt self.plot1.Run() self.y = self.y[-1] # в БД пишем решение в конце интервала self.State = sol[-1,:] def Finalize(self): self.plot1.Finalize() Задача отличается от предыдущей только составом кодом метода Run. В нём (первые 3 строчки) для решения ДУ вызывается функция odeint, а затем рассчитывается выход нашей системы (вектор y), который потом отправляется на график и в базу данных. График точно такой же, как в предыдущем случае (поэтому я его не привожу). Не забываем код для запуска задачи (и поправить имя задачи в core.py на odeint): :: model = ODE(name='odeint', TaskList=kernel.Tasks, DB=kernel.DB) model.Manager.Loop() Модель и графика разделены ~~~~~~~~~~~~~~~~~~~~~~~~~~ В этом примере (ex2.3) модель и графики выделены на разные задачи. Такую технологию можно использовать, когда у вас имеется ресурсоёмкая графика - например, трёхмерная или что-то самостоятельно написанное. Модель такая же, как в примере 2.1, и она просто разделена на две части. Задача Satellite является управляющей (isSheduler=True в создании экземпляра класса модели), её нужно запускать последней. Непрерывный объект и дискретный регулятор ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ В этом примере (ex.3) мы спроектируем цифровую систему управления для управления угловым положением космического аппарата (КА) с одной степенью свободы. Он описывается системой уравнений: .. math:: \begin{cases} \dfrac{d \omega}{dt} &= \dfrac{k \cdot M}{J},\\ \dfrac{d \varphi}{dt} &= \omega. \end{cases} Здесь :math:`\omega` - угловая скорость, :math:`M` - приложенный момент, :math:`k`- коэффициент устройства, которое приводит во вращение КА (магнитная катушка или маховик), :math:`J` - момент инерции, :math:`\varphi` - угол поворота КА. Для 1U кубсата :math:`\dfrac{k \cdot M}{J} \approx 80`, и передаточная функция будет выглядеть так: .. math:: W_{\varphi-u}(s) = \dfrac{80}{s^2}. Построим сначала корректирующее устройство в непрерывном времени, а потом дискретизируем. Буду использовать библиотеку python-control. Несколько раз пересчитывал корректирующее устройство (как это обычно получается), и вот такое получилось КУ. Код такой :: G = ct.tf([80], [1, 0, 0], inputs='u', outputs='y') T = 1/np.array([0.3, 0.0007, 30]) numR = np.polymul([T[0], 1], [T[0], 1]) denR = np.polymul([T[1], 1], [T[2], 1]) R = ct.tf(numR, denR) ct.bode_plot([G, R*G], display_margins=True, label=['$G(s)$', '$R(s)\cdot G(s)$'], dB=True, grid=True) .. image:: img/ex3_lahRGs.png ЛАХ нескорректированной системы :math:`G(s)` имеет наклон -2. Помня о том, что - переход к дискретизации будет вносить в систему неустойчивость и потянет ФЧХ вниз - желательно использовать собственную динамику системы, не тратя лишнюю энергию на ускорение делаем вывод, что вправо частоту среза переносить нецелесообразно. Проще всего разместить её вблизи :math:`\omega_c` нескорректированной системы. И даже можно немного зарезать полосу. Желаемый ЛАХ будет пересекать ось частот под наклоном -1, сопрягаться с низкочастотным участком исходной ЛАХ под -3 наклоном, а на высоких частотах пусть наклоны совпадают. Всё делается как обычно. Постоянные времени подобрал руками. Характеристики полученной системы :: closedSys = ct.feedback(R*G) S = ct.step_info(closedSys) print('Характеристики замкнутой системы:') for k in S: print(f"{k}: {S[k]}") >> Характеристики замкнутой системы: RiseTime: 1.4460521514618077 SettlingTime: 17.88236769480988 SettlingMin: 0.8554960298969676 SettlingMax: 1.405891817252921 Overshoot: 40.58918172529209 Undershoot: 0 Peak: 1.405891817252921 PeakTime: 4.037492145665642 SteadyStateValue: 1.0 Перерегулирование в 40% и время переходного порядка 17 секунд нас вполне устраивает. В библиотеке python-control есть отличный `пример дискретизации системы `_; его идея в том, что сначала мы дискретизируем и непрерывную часть системы и регулятор с частотой дискретизации большей, чем частота регулятора (в n раз), а потом просто работаем в дискретном времени. Регулятор при этом пропускает n шагов по времени, работая на своей штатной частоте (между шагами он удерживает постоянное значение на выходе, для чего мы модифицируем его методы - см. функцию sampled_data_controller). Пример этот хорош ещё и тем, что показывает способ работы и с нелинейными системами. Время шага контроллера 0.1с, дискретизируем с интервалом в n=5 раз меньше: :: Ts = .1 # sampling interval of controller dtSim = Ts/5 # time step for numerical simulation ("numerical integration") Gz = ct.c2d(G, Ts, 'zoh') Rz = ct.sample_system(ct.tf2ss(R), Ts, method='zoh', inputs='e', outputs='u') # create model of controller with a much shorter sampling time for simulation RzSim = sampled_data_controller(Rz, dtSim) GzSim = ct.c2d(G, dtSim, 'zoh') ct.bode_plot([Gz, Rz*Gz], display_margins=True, label=['$G(s)$', '$R(s)\cdot G(s)$'], dB=True, grid=True) ЛАЧХ дискретизированной системы: .. image:: img/ex3_lahRGsZ.png Видим, что система устойчива. Замкнём её, и построим переходной процесс: :: Summator = ct.summing_junction(inputs=['-y', 'r'], outputs='e') GClosedSim = ct.interconnect([RzSim, GzSim, Summator], inputs='r', outputs=['y', 'u']) # входное воздействие - скачок time = np.arange(0, 20, dtSim) step_input = np.ones_like(time) t, y = ct.input_output_response(GClosedSim, time, step_input) y, u = y # extract respones fig, axs = plt.subplots(2) axs[0].plot(t, y, '.-', label='y') axs[0].legend() axs[1].plot(t, u, '.-', label='u') axs[1].legend(); .. image:: img/ex3_Transient.png Решим ту же самую задачу, используя наш пакет. Код задачи, реализующий модель космического аппарата (в форме передаточной функции :math:`W_{\varphi-u}(s) = \dfrac{80}{s^2}`): :: import control as ct import numpy as np import libKernel import db from libGraph2D import Plot2D class Satellite(libKernel.TaskTemplate): """ модель одномерного вращательного движения КА, без трения """ def Setup(self): self.Np = 25 self.Num = [80] self.Den = [1, 0, 0] self.System = ct.tf2ss(self.Num, self.Den) self.plot1 = Plot2D(DB=self) self.plot1.Setup('time', [['y'], ['u']], chunkLen=self.Np) def Initialize(self): # задаём начальное состояние (нулевые НУ, но можно задавать из БД) self.t0 = 0. self.tSol = np.linspace(0, self.dt, self.Np) # self.State = np.zeros((len(self.Den)-1, 1)) self.State = np.zeros((self.System.A.shape[0], 1)) self.plot1.Initialize() # инициализируем график def Run(self): _, self.y, self.State = ct.forced_response(self.System, T=self.tSol, U=self.u*np.ones((1, self.tSol.shape[0])), #B.shape(1) return_x=True, squeeze=True, X0=self.State[:,-1]) self.time = self.tSol + self.t0 self.t0 += self.dt self.plot1.Run() self.y = self.y[-1] # в БД пишем решение в конце интервала def Finalize(self): self.plot1.Finalize() #% Главный цикл model = Satellite(TaskList=db.Tasks, DB=db.DataBase()) model.Manager.Loop() Видно, что эта задача включает в себя графики переходного процесса. Для того, чтобы они были детальные, используются промежуточные микрошаги (25 на 1 шаг модели dt). Регулятор реализуем по классическим формулам: .. math:: \begin{eqnarray} x[k+1] &=& A \cdot x[k] + B \cdot u[k] \\ y[k] & = & C \cdot x[k] + D \cdot u[k], \end{eqnarray} Матрицы возьмём из свойств регулятора: :: print(Rz) : sys[1]$sampled Inputs (1): ['e'] Outputs (1): ['u'] States (2): ['x[0]', 'x[1]'] A = [[ 0.0497649 0.00665116] [-0.00316722 0.99995217]] B = [[-0.03167217] [ 0.00022775]] C = [[6.86016333 0.161 ]] D = [[0.23333333]] dt = 0.1 Код задачи (tReg.py): :: import libKernel import db import numpy as np class Regulator(libKernel.TaskTemplate): """ регулятор """ def Initialize(self): # задаём начальное состояние (нулевые НУ) self.x = np.zeros((2,1)) # state self.u = 0. self.A = np.array([[ 0.0497649 , 0.00665116], [-0.00316722, 0.99995217]]) self.B = np.array([[-0.03167217], [ 0.00022775]]) self.C = np.array([[6.86016333, 0.161 ]]) self.D = np.array([[0.23333333]]) def Run(self): E = np.array([[1. - self.y]]) y = np.matmul(self.C, self.x) + np.matmul(self.D, E) self.u = float(np.squeeze(y)) self.x = np.matmul(self.A, self.x) + np.matmul(self.B, E) #% Главный цикл model = Regulator(TaskList=db.Tasks, DB=db.DataBase(), isSheduler=True) model.Manager.Loop() И база данных с описанием задач (из kernel.py): :: class DataBase(): # --- общие переменные --- t = 0. # время dt = 0.1 # шаг задач tmax = 10. # время моделирования cmd = 0 # команда всем задачам y = 0.0 # переходной процесс u = 0.0 # управление Tasks = { 'Satellite': {'Keys':'t,cmd,dt,y,u'}, 'Regulator': {'Keys':'t,cmd,dt,y,u'}, } Переходной процесс такой же, как и полученный с помощью python-control: .. image:: img/ex3_Results.png Видно только задержку на один такт симуляции, в остальном же переходной процесс полностью совпадает.