Использование¶
Установка¶
В этой версии нужно иметь следующие файлы (в каталоге, куда настроены пути импорта):
libKernel.py
libComm.py
libGraph2D.py
Зависимости¶
matplotlib
numpy
pyquaternion
Для запуска примеров:
scipy
control
cobs
pyserial
Установите их своим пакетным менеджером (mamba, conda или pip):
(.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
Примеры¶
Рассмотрим несколько вариантов моделей и сравним их с аналогичными моделями, выполненными в SimInTech и с теоретически рассчитанными переходными процессами.
Передаточная функция¶
Смоделируем переходной процесс в передаточной функции (ПФ) при нулевых начальных условиях от входного воздействия типа единичный скачок. Для понимания кода нужно иметь опыт работы с библиотекой pycontrol.
Используем ПФ вида
Обозначим входное воздействие переменной 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()
Рассчитаем переходной процесс двумя способами:
Разобьём наш интервал времени на шаги равные dt, и будем считать решение только в этих точках.
Каждый из шагов разобьём на микрошаги, чтобы получить более детальный переходной процесс.
Способ с большими шагами¶
Код задачи для первого способа будет выглядеть следующим образом (пример 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 график финализируется - при окончании расчёта задачи график дорисовывает все точки, которые он успел получить.
В итоге откроется такой график:
Микрошаги¶
Это способ с промежуточными шагами (пример 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, но с микрошагами:
Совместив оба графика, видим, что точки решения совпадают в моменты времени, кратные dt:
Значение выходной переменной y и при большом шаге, и при разбивке на микрошаги получилось одинаковым - потому что библиотека control для систем с непрерывным временем использует аналитический расчёт ПП с использованием матричной экспоненты.
Дифференциальное уравнение¶
В случае, когда требуется решить дифференциальное уравнение общего вида, удобно применить функцию odeint библиотеки scipy.
Наша передаточная функция в форме Коши (или в форме пространства состояний) будет иметь вид
где
Зададим правую часть ДУ:
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)
Код задачи выглядит следующим образом (пример 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) мы спроектируем цифровую систему управления для управления угловым положением космического аппарата (КА) с одной степенью свободы. Он описывается системой уравнений:
Здесь \(\omega\) - угловая скорость, \(M\) - приложенный момент, \(k\)- коэффициент устройства, которое приводит во вращение КА (магнитная катушка или маховик), \(J\) - момент инерции, \(\varphi\) - угол поворота КА.
Для 1U кубсата \(\dfrac{k \cdot M}{J} \approx 80\), и передаточная функция будет выглядеть так:
Построим сначала корректирующее устройство в непрерывном времени, а потом дискретизируем. Буду использовать библиотеку 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)
ЛАХ нескорректированной системы \(G(s)\) имеет наклон -2. Помня о том, что
переход к дискретизации будет вносить в систему неустойчивость и потянет ФЧХ вниз
желательно использовать собственную динамику системы, не тратя лишнюю энергию на ускорение
делаем вывод, что вправо частоту среза переносить нецелесообразно. Проще всего разместить её вблизи \(\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)
ЛАЧХ дискретизированной системы:
Видим, что система устойчива. Замкнём её, и построим переходной процесс:
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();
Решим ту же самую задачу, используя наш пакет.
Код задачи, реализующий модель космического аппарата (в форме передаточной функции \(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).
Регулятор реализуем по классическим формулам:
Матрицы возьмём из свойств регулятора:
print(Rz)
<StateSpace>: 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:
Видно только задержку на один такт симуляции, в остальном же переходной процесс полностью совпадает.