Постановка задач

Решётчатые функции определены в дискретные моменты времени, которые, в общем случае, могут находиться на разном расстоянии друг от друга. Пусть h_p — шаг решётки, тогда, в общем случае, $h_p\neq const$ . Для нас наиболее интересна ситуация когда h_p = const . Это соответствует характеру работы аналого-цифрового преобразователя (АЦП) управляющего цифрового устройства, когда значения технологического параметра преобразуются в цифровой формат с некоторой дискретностью. Последовательность таких значений описывает непрерывную функцию технологического параметра в решётчатом виде.

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

Выделим следующие задачи:

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

Выбор метрики

На практике широко используются относительная и приведённая погрешности.

Относительная погрешность чаще всего рассчитывается в процентах:

$\delta Y=\frac{Y_{изм}-Y_{ист}}{Y_{ист}}\cdot 100\%$

где $Y_{изм}$ — измеренное значение параметра в текущий момент времени, $Y_{ист}$ — истинное значение параметра в этот момент времени.

Для подсчитанных значений относительной погрешности можно рассчитать среднюю относительную погрешность на данных измерениях:

$\delta Y_{mean}=\delta Y_m=\frac{100\%}{n}\sum^n_{i=1}\frac{Y_{изм_i}-Y_{ист_i}}{Y_{ист_i}}$

Назовем среднюю относительную погрешность метрикой MLE — Mean reLative Error.

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

Приведённая погрешность также чаще всего вычисляется в процентах:

$\gamma Y=\frac{Y_{изм}-Y_{ист}}{Y_{норм}}\cdot 100\%$

где $Y_{изм}$ — измеренное значение параметра в текущий момент времени, $Y_{норм}$ — нормирующее значение, равное диапазону изменения измеряемой величины или, иными словами, шкале измерительного прибора.

Средняя приведённая погрешность:

$\gamma Y_{mean}=\gamma Y_m=\frac{100\%}{n}\sum^n_{i=1}\frac{Y_{изм_i}-Y_{ист_i}}{Y_{норм}}$

Назовем среднюю приведённую погрешность $\gamma Y_m$ метрикой MRE — Mean Reduced Error.

MRE, как и MLE, имеет ясный физический смысл — она показывает насколько в среднем в процентах отклонилась прогнозная / измеренная величина от диапазона изменения / шкалы прибора истинной величины. Величина $Y_{норм}=Y_{ист_{max}}-Y_{ист_{min}}$ всегда больше нуля. Поэтому MRE лишена недостатков MLE и может быть смело использована для оценки качества прогноза.

Скрытый текст

def MRE(y_true, y_pred, arr=False):
    """
    Средняя приведенная к диапазону изменения переменной ошибка

    y_true - истинное значение
    y_pred - прогнозное значение
    arr   - массив значений абсолютной ошибки, в долях единицы
            к истинному значению:
     True - вернуть
    False - не возвращать
    
    По умолчанию возвращает среднюю абсолютную процентную ошибку, %
    """

    y_t = list(y_true)
    y_p = list(y_pred)

    interval = max(y_t) - min(y_t)

    mre = []

    for i in range(len(y_t)):
        mre.append(np.abs((y_t[i] - y_p[i]) / interval))

    if arr:
        return mre
    else:
        if np.isinf(np.mean(mre)):
            return BIGNUM
        else:
            return np.mean(mre) * 100

Моделирование типовых динамических звеньев

# Импорт библиотек
import math
import sys
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import control as ct

BIGNUM = sys.maxsize

s = ct.tf('s') # Оператор передаточной функции в преобразовании Лапласа

mt = 100 # Время моделирования, с
hp = 0.1 # Шаг решетчатой функции, с
mtl = [i * hp for i in np.arange(int(math.ceil((mt + hp)/hp)))] # ВременнАя шкала моделирования, с

Апериодическое звено 1-го порядка

$W(s) = \frac{k}{Ts + 1} = \frac{1}{15s+1}$ $T\,\frac{dy(t)}{dt} + y(t) = k\,x(t)$ $\boxed{15\,\frac{dy(t)}{dt} + y(t) = x(t)},\,k=1,\,T=15$

k = 1
T = 15

A1 = (k) / (T * s + 1)

t, y_A1 = ct.step_response(A1, T=mtl)

Апериодическое звено 2-го порядка

$W(s) = \frac{k}{T_{2}^2s^2 + T_{1}s + 1} = \frac{1}{100s^2+20s+1},\,\,T_1\,\geqslant\,2T_2$ $T_2^2\,\frac{d^2y(t)}{dt} + T_1\,\frac{dy(t)}{dt} + y(t) = k\,x(t),$ $\boxed{100\,\frac{d^2y(t)}{dt} + 20\,\frac{dy(t)}{dt} + y(t) = x(t)},\,k=1,\,T_1=20,\,T_2=10$

k = 1
T1 = 20
T2 = 10

A2 = (k) / (T2**2 * s**2 + T1 * s + 1)

t, y_A2 = ct.step_response(A2, T=mtl)

Колебательное звено

$W(s) =\frac{k}{T^2s^2 + 2\xi Ts + 1} = \frac{1}{16s^2+2s+1},\,\,0<\xi<1$ $T^2\,\frac{d^2y(t)}{dt} + 2\,\xi\,T\,\frac{dy(t)}{dt} + y(t) = k\,x(t),$ $\boxed{16\,\frac{d^2y(t)}{dt} + 2\,{\cdot}\,0,25\,{\cdot}\,4\frac{dy(t)}{dt} + y(t) = x(t)},\,k=1,\,T=4,\,\xi=0,25$

k = 1
T = 4
ksi = 0.25

Ksi = (k) / (T**2 * s**2 + 2 * ksi * T * s + 1)

t, y_Ksi = ct.step_response(Ksi, T=mtl)

Консервативное звено

$W(s) = \frac{k}{T^2s^2 + 1} = \frac{0.75}{16s^2+1}$ $T^2\,\frac{d^2y(t)}{dt} + y(t) = k\,x(t),$ $\boxed{4\,\frac{d^2y}{dt} + y = x(t)},\,k=0.75,\,T=4$

k = 0.75
T = 4

Cons = (k) / (T**2 * s**2 + 1)

t, y_Cons = ct.step_response(Cons, T=mtl)

Аддитивный белый гауссовский шум

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

Ступенчатое гауссовское возмущение

В качестве возмущающего будем применять стахостическое (гауссовское) ступенчатое воздействие. Назовём его функцией Хевисайда-Гаусса.

Моделирование линейной регрессией типовых динамических звеньев

Динамические свойства звеньев будем описывать соответствующими дифференциальными уравнениями (ДУ). Затем будем переходить к формам ДУ для решётчатых функций и строить на их основе итерируемые схемы прогнозирования очередного значения функции.

Для нахождения дифференциалов решётчатых функций широко применяются численные методы.

Дифференциалы (производные) решётчатой функции будем обозначать, как $\hat y', \,\hat y''$ .

Используя формулу Лагранжа для случая равных промежутков формула численного дифференцирования производной первого порядка с четвёртым порядком погрешности примет вид:

$\hat y'_i = \frac{1}{12\,h_P}\left(3y_{i-4} - 16\,y_{i-3} + 36\,y_{i-2} - 48\,y_{i-1} + 25\,y_i\right)$

Производная второго порядка с третьим порядком погрешности примет вид:

$\hat y''_i = \frac{1}{12\,h^2_P}\left(11y_{i-4} - 56\,y_{i-3} + 114\,y_{i-2} - 104\,y_{i-1} + 35\,y_i\right)$

Апериодическое звено 1-го порядка

А-звено 1-го порядка (А1-звено) описывается следующим дифференциальным уравнением:

$T\,\frac{dy(t)}{dt} + y(t) = k\,x(t),\; T\,\frac{dy}{dt} + y = k\,x$

Для решётчатой функции данное уравнение примет вид:

$T\hat y' + y = k\,x$

Применяя формулы Тейлора, для -го значения решётчатой функции дифференциальное уравнение звена примет вид:

$T{\cdot}\hat y'_i + y_i = k{\cdot}x_i$

Тогда:

$\frac{T}{12\,h_P}{\cdot}(3\,y_{i-4} - 16\,y_{i-3} + 36\,y_{i-2} - 48\,y_{i-1} + 25\,y_i) + y_i = k{\cdot}x_i$ $3\,y_{i-4} - 16\,y_{i-3} + 36\,y_{i-2} - 48\,y_{i-1} + 25\,y_i + \frac{12\,h_P}{T}{\cdot}y_i = \frac{12\,h_P}{T}{\cdot}k{\cdot}x_i$ $y_i = \frac{12\,k\,h_P}{25\,T + 12\,h_P}{\cdot}x_i\,-\,\frac{3\,y_{i-4} - 16\,y_{i-3} + 36\,y_{i-2} - 48\,y_{i-1}}{25\,+\,\frac{12\,h_P}{T}}$ $y_i = C_1{\cdot}x_i\,+\,C_2{\cdot}y_{i-1}\,-\,C_3{\cdot}y_{i-2}\,+\, C_4{\cdot}y_{i-3}\,-\,C_5{\cdot}y_{i-4}$ $C_1 = \frac{12\,k\,h_P}{25\,T + 12\,h_P},\,\, C_2 = \frac{48}{25\,+\,\frac{12\,h_P}{T}},\,\, C_3 = \frac{36}{25\,+\,\frac{12\,h_P}{T}}$ $C_4 = \frac{16}{25\,+\,\frac{12\,h_P}{T}},\,\, C_5 = \frac{3}{25\,+\,\frac{12\,h_P}{T}}$

hp = 0.1
k = 1
T = 15

x_t = [next(pg) for i in mtl]

t, y_A1 = ct.forced_response(A1, T=mtl, U=x_t)

C1 = 12 * k * hp / (25 * T + 12 * hp)
C2 = 48 / (25 + 12 * hp / T)
C3 = 36 / (25 + 12 * hp / T)
C4 = 16 / (25 + 12 * hp / T)
C5 = 3 / (25 + 12 * hp / T)

y_A1_lr = [0, 0, 0, 0]  # Нулевые начальные условия

for i in range(0, len(x_t)):
    y_A1_lr.append(
        C1 * x_t[i] + C2 * y_A1_lr[-1] - C3 * y_A1_lr[-2] + C4 * y_A1_lr[-3] - C5 * y_A1_lr[-4]
    )

y_A1_lr = y_A1_lr[4:]

Апериодическое звено 2-го порядка

А-звено 2-го порядка описывается (А2-звено) следующим дифференциальным уравнением:

$T_2^2\,\frac{d^2y(t)}{dt} + T_1\,\frac{dy(t)}{dt} + y(t) = k\,x(t)$ $T_2^2\,\frac{d^2y}{dt} + T_1\,\frac{dy}{dt} + y = k\,x, \,\,T_1 \geqslant 2\, T_2 \; \left(T_2 \leqslant \frac{T_1}{2}\right)$

Для решётчатой функции данное уравнение примет вид:

$T_2^2\hat y'' + T_1\hat y' + y = k\,x$

Для решётчатой функции данное уравнение примет вид:

$T_2^2\hat y'' + T_1\hat y' + y = k\,x$

Тогда:

$\frac{T_2^2}{12\,h^2_P}{\cdot}(11\,y_{i-4} - 56\,y_{i-3} + 114\,y_{i-2} - 104\,y_{i-1} + 35\,y_i) +$ $+\frac{T_1}{12\,h_P}{\cdot}(3\,y_{i-4} - 16\,y_{i-3} + 36\,y_{i-2} - 48\,y_{i-1} + 25\,y_i) + y_i = k{\cdot}x_i$ $\frac{11\,T_2^2}{12h^2_P}{\cdot}y_{i-4} - \frac{56\,T_2^2}{12h^2_P}{\cdot}y_{i-3} + \frac{114\,T_2^2}{12h^2_P}{\cdot}y_{i-2} - \frac{104\,T_2^2}{12h^2_P}{\cdot}y_{i-1} + \frac{35\,T_2^2}{12h^2_P}{\cdot}y_i +$ $+ \frac{3\,T_1}{12h_P}{\cdot}y_{i-4} - \frac{16\,T_1}{12h_P}{\cdot}y_{i-3} + \frac{36\,T_1}{12h_P}{\cdot}y_{i-2} - \frac{48\,T_1}{12h_P}{\cdot}y_{i-1} + \frac{25\,T_1}{12h_P}{\cdot}y_i + y_i = k\,x_i$ $\left(\frac{11\,T_2^2}{12\,h^2_P} +\frac{3\,T_1}{12\,h_P}\right){\cdot}y_{i-4} - \left(\frac{56\,T_2^2}{12\,h^2_P} + \frac{16\,T_1}{12\,h_P}\right){\cdot}y_{i-3} + \left(\frac{114\,T_2^2}{12\,h^2_P} + \frac{36\,T_1}{12\,h_P}\right){\cdot}y_{i-2} -$ $- \left(\frac{104\,T_2^2}{12\,h^2_P} + \frac{48\,T_1}{12\,h_P}\right){\cdot}y_{i-1} + \left(\frac{35\,T_2^2}{12h^2_P} +\frac{25\,T_1}{12\,h_P} + 1\right){\cdot}y_i = k{\cdot}x_i$ $y_i = C_1{\cdot}x_i\,+\,C_2{\cdot}y_{i-1}\,-\,C_3{\cdot}y_{i-2}\,+\, C_4{\cdot}y_{i-3}\,-\,C_5{\cdot}y_{i-4}$ $C_1 = \frac{k}{\frac{35\,T_2^2}{12\,h^2_P} +\frac{25\,T_1}{12\,h_P} + 1},\,\, C_2 = \frac{\frac{104\,T_2^2}{12\,h^2_P} + \frac{48\,T_1}{12\,h_P}}{\frac{35\,T_2^2}{12\,h^2_P} +\frac{25\,T_1}{12\,h_P} + 1},\,\, C_3 = \frac{\frac{114\,T_2^2}{12\,h^2_P} + \frac{36\,T_1}{12\,h_P}}{\frac{35\,T_2^2}{12\,h^2_P} +\frac{25\,T_1}{12\,h_P} + 1}$ $C_4 = \frac{\frac{56\,T_2^2}{12\,h^2_P} + \frac{16\,T_1}{12\,h_P}}{\frac{35\,T_2^2}{12\,h^2_P} +\frac{25\,T_1}{12\,h_P} + 1},\,\, C_5 = \frac{\frac{11\,T_2^2}{12\,h^2_P} +\frac{3\,T_1}{12\,h_P}}{\frac{35\,T_2^2}{12\,h^2_P} +\frac{25\,T_1}{12\,h_P} + 1}$

k = 1
T1 = 20
T2 = 10

x_t = [next(pg) for i in mtl]

t, y_A2 = ct.forced_response(A2, T=mtl, U=x_t)

_d = (35 * T2**2) / (12 * hp**2) + (25 * T1) / (12 * hp) + 1

C1 = k / _d
C2 = ((104 * T2**2) / (12 * hp**2) + (48 * T1) / (12 * hp)) / _d
C3 = ((114 * T2**2) / (12 * hp**2) + (36 * T1) / (12 * hp)) / _d
C4 = ((56 * T2**2) / (12 * hp**2) + (16 * T1) / (12 * hp)) / _d
C5 = ((11 * T2**2) / (12 * hp**2) + (3 * T1) / (12 * hp)) / _d

y_A2_lr = [0, 0, 0, 0]  # Нулевые начальные условия

for i in range(0, len(x_t)):
    y_A2_lr.append(
        C1 * x_t[i] + C2 * y_A2_lr[-1] - C3 * y_A2_lr[-2] + C4 * y_A2_lr[-3] - C5 * y_A2_lr[-4]
    )

y_A2_lr = y_A2_lr[4:]

Колебательное звено

Колебательное звено ( $\xi$ -звено) описывается следующим дифференциальным уравнением:

$T^2\,\frac{d^2y(t)}{dt} + 2\,\xi\,T\,\frac{dy(t)}{dt} + y(t) = k\,x(t)$ $T^2\,\frac{d^2y}{dt} + 2\,\xi\,T\,\frac{dy}{dt} + y = k\,x,\,0<\xi<1$

Для решётчатой функции данное уравнение примет вид:

$T^2\,\hat y'' + 2\,\xi\,T\,\hat y' + y = k\,x$

Применяя формулы Тейлора, для -го значения решётчатой функции дифференциальное уравнение звена примет вид:

$T^2\,\hat y_i'' + 2\,\xi\,T\,\hat y_i' + y_i = k\,x$

Тогда:

$\frac{T^2}{12\,h^2_P}{\cdot}(11\,y_{i-4} - 56\,y_{i-3} + 114\,y_{i-2} - 104\,y_{i-1} + 35\,y_i) +$ $+\frac{2\,\xi\,T}{12\,h_P}{\cdot}(3\,y_{i-4} - 16\,y_{i-3} + 36\,y_{i-2} - 48\,y_{i-1} + 25\,y_i) + y_i = k{\cdot}x_i$ $\frac{11\,T^2}{12\,h^2_P}{\cdot}y_{i-4} - \frac{56\,T^2}{12\,h^2_P}{\cdot}y_{i-3} + \frac{114\,T^2}{12\,h^2_P}{\cdot}y_{i-2} - \frac{104\,T^2}{12\,h^2_P}{\cdot}y_{i-1} + \frac{35\,T^2}{12\,h^2_P}{\cdot}y_i +$ $+ \frac{6\,\xi\,T}{12\,h_P}{\cdot}y_{i-4} - \frac{32\,\xi\,T}{12\,h_P}{\cdot}y_{i-3} + \frac{72\,\xi\,T}{12\,h_P}{\cdot}y_{i-2} - \frac{96\,\xi\,T}{12\,h_P}{\cdot}y_{i-1} + \frac{50\,\xi\,T}{12\,h_P}{\cdot}y_i + y_i = k{\cdot}x_i$ $\left(\frac{11\,T^2}{12\,h^2_P} +\frac{6\,\xi\,T}{12\,h_P}\right){\cdot}y_{i-4} - \left(\frac{56\,T^2}{12\,h^2_P} + \frac{32\,\xi\,T}{12\,h_P}\right){\cdot}y_{i-3} + \left(\frac{114\,T^2}{12\,h^2_P} + \frac{72\,\xi\,T}{12\,h_P}\right){\cdot}y_{i-2} -$ $- \left(\frac{104\,T^2}{12\,h^2_P} + \frac{96\,\xi\,T}{12\,h_P}\right){\cdot}y_{i-1} + \left(\frac{35\,T^2}{12\,h^2_P} +\frac{50\,\xi\,T}{12\,h_P} + 1\right){\cdot}y_i = k{\cdot}x_i$ $y_i = C_1{\cdot}x_i\,+\,C_2{\cdot}y_{i-1}\,-\,C_3{\cdot}y_{i-2}\,+\, C_4{\cdot}y_{i-3}\,-\,C_5{\cdot}y_{i-4}$ $C_1 = \frac{k}{\frac{35\,T^2}{12\,h^2_P} +\frac{50\,\xi\,T}{12\,h_P} + 1},\,\, C_2 = \frac{\frac{104\,T^2}{12\,h^2_P} + \frac{96\,\xi\,T}{12\,h_P}}{\frac{35\,T^2}{12\,h^2_P} +\frac{50\,\xi\,T}{12\,h_P} + 1},\,\, C_3 = \frac{\frac{114\,T^2}{12\,h^2_P} + \frac{72\,\xi\,T}{12\,h_P}}{\frac{35\,T^2}{12\,h^2_P} +\frac{50\,\xi\,T}{12\,h_P} + 1}$ $C_4 = \frac{\frac{56\,T^2}{12\,h^2_P} + \frac{32\,\xi\,T}{12\,h_P}}{\frac{35\,T^2}{12\,h^2_P} +\frac{50\,\xi\,T}{12\,h_P} + 1},\,\, C_5 = \frac{\frac{11\,T^2}{12\,h^2_P} +\frac{6\,\xi\,T}{12\,h_P}}{\frac{35\,T^2}{12\,h^2_P} +\frac{50\,\xi\,T}{12\,h_P} + 1}$

k = 1
T = 4
ksi = 0.25

x_t = [next(pg) for i in mtl]

t, y_Ksi = ct.forced_response(Ksi, T=mtl, U=x_t)

_d = (35 * T**2) / (12 * hp**2) + (50 * ksi * T) / (12 * hp) + 1 

C1 = k / _d
C2 = ((104 * T**2) / (12 * hp**2) + (96 * ksi * T) / (12 * hp)) / _d
C3 = ((114 * T**2) / (12 * hp**2) + (72 * ksi * T) / (12 * hp)) / _d
C4 = ((56 * T**2) / (12 * hp**2) + (32 * ksi * T) / (12 * hp)) / _d
C5 = ((11 * T**2) / (12 * hp**2) + (6 * ksi * T) / (12 * hp)) / _d

y_Ksi_lr = [0, 0, 0, 0]  # Нулевые начальные условия

for i in range(0, len(x_t)):
    y_Ksi_lr.append(
        C1 * x_t[i] + C2 * y_Ksi_lr[-1] - C3 * y_Ksi_lr[-2] + C4 * y_Ksi_lr[-3] - C5 * y_Ksi_lr[-4]
    )

y_Ksi_lr = y_Ksi_lr[4:]

Консервативное звено

Консервативное звено (-звено) описывается следующим дифференциальным уравнением:

$T^2\,\frac{d^2y(t)}{dt} + y(t) = k\,x(t)$ $T^2\,\frac{d^2y}{dt} + y = k\,x$

Для решётчатой функции данное уравнение примет вид:

$T^2\hat y'' + y = k\,x$

Применяя формулы Тейлора, для -го значения решётчатой функции дифференциальное уравнение звена примет вид:

$T^2{\cdot}\hat y''_i + y_i = k{\cdot}x_i$

Тогда:

$\frac{T^2}{12\,h^2_P}{\cdot}(11\,y_{i-4} - 56\,y_{i-3} + 114\,y_{i-2} - 104\,y_{i-1} + 35\,y_i) + y_i = k{\cdot}x_i$ $11\,y_{i-4} - 56\,y_{i-3} + 114\,y_{i-2} - 104\,y_{i-1} + 35\,y_i + \frac{12\,h^2_P}{T^2}{\cdot}y_i = \frac{12\,h^2_P}{T^2}{\cdot}k{\cdot}x_i$ $y_i = \frac{12\,k\,h^2_P}{35\,T^2 + 12\,h^2_P}{\cdot}x_i\,-\,\frac{11\,y_{i-4} - 56\,y_{i-3} + 114\,y_{i-2} - 104\,y_{i-1}}{35\,+\,\frac{12\,h^2_P}{T^2}}$ $y_i = C_{1}{\cdot}x_i\,+\,C_{2}{\cdot}y_{i-1}\,-\,C_{3}{\cdot}y_{i-2}\,+\, C_{4}{\cdot}y_{i-3}\,-\,C_{5}{\cdot}y_{i-4}$ $C_{1} = \frac{12\,k\,h^2_P}{35\,T^2 + 12\,h^2_P},\,\, C_{2} = \frac{104}{35\,+\,\frac{12\,h^2_P}{T^2}},\,\, C_{3} = \frac{114}{35\,+\,\frac{12\,h^2_P}{T^2}}$ $C_{4} = \frac{56}{35\,+\,\frac{12\,h^2_P}{T^2}},\,\, C_{5} = \frac{11}{35\,+\,\frac{12\,h^2_P}{T^2}}$

k = 0.75
T = 4

x_t = [next(pg) for i in mtl]

t, y_Cons = ct.forced_response(Cons, T=mtl, U=x_t)

C1 = 12 * k * hp**2 / (35 * T**2 + 12 * hp**2)
C2 = 104 / (35 + 12 * hp**2 / T**2)
C3 = 114 / (35 + 12 * hp**2 / T**2)
C4 = 56 / (35 + 12 * hp**2 / T**2)
C5 = 11 / (35 + 12 * hp**2 / T**2)

y_Cons_lr = [0, 0, 0, 0]  # Нулевые начальные условия

for i in range(0, len(x_t)):
    y_Cons_lr.append(
        C1 * x_t[i] + C2 * y_Cons_lr[-1] - C3 * y_Cons_lr[-2] + C4 * y_Cons_lr[-3] - C5 * y_Cons_lr[-4]
    )

y_Cons_lr = y_Cons_lr[4:]

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

# A1-звено
MRE(y_A1, y_A1_lr)

0.037048140997006936

# A2-звено
MRE(y_A2, y_A2_lr)

0.04460863218205424

#Ksi-звено
MRE(y_Ksi, y_Ksi_lr)

0.013608750809221084

#Cons-звено
MRE(y_Cons, y_Cons_lr)

0.026759652627513193

Поскольку на практике чаще встречаются технологические параметры, описываемые ДУ второго и более высоких порядков и поскольку колебательное и консервативное звенья являются частными случаями А2-звена, далее будем исследовать А2-звено.

Определение коэффициентов линейной регрессии по временным рядам решётчатой функции

Поиск параметров $C_1,\,... ,\,C_5$ аналитическим способом если и возможен, то крайне затруднителен. В конечном итоге он сводится к подбору таких значений $C_1,\,... ,\,C_5$ , при которых MRE минимальна. Иными словами:

$MRE = f(C_1, C_2, C_3, C_4, C_5) = f_{MRE}$

и наша задача — $arg\,min f_{MRE}$ .

Параметры $C_1,\,... ,\,C_5$ рассчитываются как комбинации параметров звена, в случае -звена — $k,\,T_1,\,T_2$ . Тогда:

$f_{MRE} = f(C_1,C_2, C_3, C_4, C_5) = f(k,\,T_1,\,T_2)$

Задача $arg\,min f_{MRE}(k,\,T_1, T_2)$ — задача поиска глобального минимума функции трёх переменных.

Для поиска глобального минимума функции нескольких переменных воспользуемся, например, функцией minimize() библиотеки scipy.

Скрытый текст

from scipy.optimize import minimize

def f_A2(n, *args):
    '''
    Функция проверки коэффициентов A2-звена
    
    n - параметры звена
    args[0] - x - входные данные (возмущение)
    args[1] - y - выходные данные (реакция)
              hp - шаг решётчатой функции
              p - глубина истории, n-1 для n-точечного шаблона
    '''

    x = args[0].copy()
    y = args[1].copy()
    hp = 0.1
    p = 4

    k = n[0]
    T1 = n[1]
    T2 = n[2]

    _d = (35 * T2**2) / (12 * hp**2) + (25 * T1) / (12 * hp) + 1

    C1 = k / _d
    C2 = ((104 * T2**2) / (12 * hp**2) + (48 * T1) / (12 * hp)) / _d
    C3 = ((114 * T2**2) / (12 * hp**2) + (36 * T1) / (12 * hp)) / _d
    C4 = ((56 * T2**2) / (12 * hp**2) + (16 * T1) / (12 * hp)) / _d
    C5 = ((11 * T2**2) / (12 * hp**2) + (3 * T1) / (12 * hp)) / _d
    
    x_t = x
    y_t = args[2].copy()  # Начальные условия

    for i in range(0, len(x_t)):  # Учет нулевых начальных значений -- y(x_t) = 0
        y_t.append(C1 * x_t[i]
                 + C2 * y_t[-1]
                 - C3 * y_t[-2]
                 + C4 * y_t[-3]
                 - C5 * y_t[-4]
                  )
            
    y_t = y_t[p:]

    return MRE(y, y_t)

x_t = [next(pg) for i in mtl]
t, y_A2 = ct.forced_response(A2, T=mtl, U=x_t)

res = minimize(f_A2,
               [0.01, 0.01, 0.01],
               args=(x_t, y_A2, [0, 0, 0, 0]),
               method="Powell",
               options={'xtol': 1e-8, 'disp': True},
               bounds=[(0.01, 5), (0.01, 30), (0.01, 20)])

print('\nMRE:', res.fun)
print('k:', res.x[0])
print('T1:', res.x[1])
print('T2:', res.x[2])

#--------------------------------------
# Optimization terminated successfully.
#          Current function value: 0.128574
#          Iterations: 13
#          Function evaluations: 1505
# 
# MRE: 0.12857363263770344
# k: 0.9997163477435754
# T1: 19.998935106980777
# T2: 9.9991512679897

Проверим работоспособность подхода на ненулевых исходных данных.

j = 100

res = minimize(f_A2,
               [0.01, 0.01, 0.01],
               args=(x_t[j:], y_A2[j:], [y_A2[j-4], y_A2[j-3], y_A2[j-2], y_A2[j-1]]),
               method="Powell",
               options={'xtol': 1e-8, 'disp': True},
               bounds=[(0.01, 5), (0.01, 30), (0.01, 20)])

print('\nMRE:', res.fun)
print('k:', res.x[0])
print('T1:', res.x[1])
print('T2:', res.x[2])

#--------------------------------------
# Optimization terminated successfully.
#          Current function value: 0.000829
#          Iterations: 9
#          Function evaluations: 999
# 
# MRE: 0.0008285555371525098
# k: 0.9999871396033223
# T1: 19.99983473603636
# T2: 9.99994226212653

Полученный результат практически идеален.

Зашумлённые данные

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

x_t = [next(pg) for i in mtl]
t, y_A2 = ct.forced_response(A2, T=mtl, U=x_t)

ln = max(abs(max(y_A2)), abs(min(y_A2))) * 0.1

y_A2_g = [y_A2[i]+gauss()*ln for i in range(len(mtl))]

j = 100

res = minimize(f_A2,
               [0.01, 0.01, 0.01],
               args=(x_t[j:], y_A2_g[j:], [y_A2_g[j-4], y_A2_g[j-3], y_A2_g[j-2], y_A2_g[j-1]]),
               method="Powell",
               options={'xtol': 1e-8, 'disp': True},
               bounds=[(0.01, 5), (0.01, 30), (0.01, 20)])

print('\nMRE:', res.fun)
print('k:', res.x[0])
print('T1:', res.x[1])
print('T2:', res.x[2])

#--------------------------------------
# Optimization terminated successfully.
#          Current function value: 2.530822
#          Iterations: 6
#          Function evaluations: 676
# 
# MRE: 2.5308224410672238
# k: 1.014091723635805
# T1: 20.44728805919973
# T2: 10.622387452988548

Полученные результаты будут ухудшаться при росте амплитуды аддитивного гауссовского шума. Это связано с тем, что начальные условия -- четыре значения функции до момента начала прогнозирования -- не совпадают с истинными.

Обозначим $f^0_{-4},\,f^0_{-3},\,f^0_{-2},\,f^0_{-1}$ как значения начальных условий. Тогда задача $arg\,min f_{MRE}$ примет вид:

$arg\,min f_{MRE}(k,\,T_1,\,T_2,\,f^0_{-4},\,f^0_{-3},\,f^0_{-2},\,f^0_{-1})$

Скрытый текст

def f_A2_g(n, *args):
    '''
    Функция проверки коэффициентов A2-звена
    
    n - параметры звена, начальные условия
    args[0] - x - входные данные (возмущение)
    args[1] - y - выходные данные (реакция)
              hp - шаг решётчатой функции
              p - глубина истории, n-1 для n-точечного шаблона
    '''

    x = args[0].copy()
    y = args[1].copy()
    hp = 0.1
    p = 4

    k = n[0]
    T1 = n[1]
    T2 = n[2]

    _d = (35 * T2**2) / (12 * hp**2) + (25 * T1) / (12 * hp) + 1

    C1 = k / _d
    C2 = ((104 * T2**2) / (12 * hp**2) + (48 * T1) / (12 * hp)) / _d
    C3 = ((114 * T2**2) / (12 * hp**2) + (36 * T1) / (12 * hp)) / _d
    C4 = ((56 * T2**2) / (12 * hp**2) + (16 * T1) / (12 * hp)) / _d
    C5 = ((11 * T2**2) / (12 * hp**2) + (3 * T1) / (12 * hp)) / _d
    
    x_t = x
    y_t = [n[3], n[4], n[5], n[6]]  # Начальные условия

    for i in range(0, len(x_t)):  # Учет нулевых начальных значений -- y(x_t) = 0
        y_t.append(C1 * x_t[i]
                 + C2 * y_t[-1]
                 - C3 * y_t[-2]
                 + C4 * y_t[-3]
                 - C5 * y_t[-4]
                  )
            
    y_t = y_t[p:]

    return MRE(y, y_t)

x_t = [next(pg) for i in mtl]
t, y_A2 = ct.forced_response(A2, T=mtl, U=x_t)

ln = max(abs(max(y_A2)), abs(min(y_A2))) * 0.1

y_A2_g = [y_A2[i]+gauss()*ln for i in range(len(mtl))]

j = 100

max_v = max(y_A2_g[j-4],
            y_A2_g[j-3],
            y_A2_g[j-2],
            y_A2_g[j-1])

min_v = min(y_A2_g[j-4],
            y_A2_g[j-3],
            y_A2_g[j-2],
            y_A2_g[j-1])

mean_v = (min_v + max_v) / 2

res = minimize(f_A2_g,
               [0.01, 0.01, 0.01, mean_v, mean_v, mean_v, mean_v],
               args=(x_t[j:], y_A2_g[j:]),
               method="Powell",
               options={'xtol': 1e-12, 'disp': True},
               bounds=[(0.01, 5), (0.01, 30), (0.01, 20), (min_v, max_v), (min_v, max_v), (min_v, max_v), (min_v, max_v)])

print('\nMRE:', res.fun)
print('k:', res.x[0])
print('T1:', res.x[1])
print('T2:', res.x[2])

#--------------------------------------
# Optimization terminated successfully.
#          Current function value: 2.318385
#          Iterations: 8
#          Function evaluations: 2286
# 
# MRE: 2.318384907285262
# k: 0.9939295663300647
# T1: 19.81697562916168
# T2: 9.959654243923332

x_t = [next(pg) for i in mtl]
t, y_A2 = ct.forced_response(A2, T=mtl, U=x_t)

ln = max(abs(max(y_A2)), abs(min(y_A2))) * 0.5

y_A2_g = [y_A2[i]+gauss()*ln for i in range(len(mtl))]

j = 100

max_v = max(y_A2_g[j-4],
            y_A2_g[j-3],
            y_A2_g[j-2],
            y_A2_g[j-1])

min_v = min(y_A2_g[j-4],
            y_A2_g[j-3],
            y_A2_g[j-2],
            y_A2_g[j-1])

mean_v = (min_v + max_v) / 2

res = minimize(f_A2_g,
               [0.01, 0.01, 0.01, mean_v, mean_v, mean_v, mean_v],
               args=(x_t[j:], y_A2_g[j:]),
               method="Powell",
               options={'xtol': 1e-12, 'disp': True},
               bounds=[(0.01, 5), (0.01, 30), (0.01, 20), (min_v, max_v), (min_v, max_v), (min_v, max_v), (min_v, max_v)])

print('\nMRE:', res.fun)
print('k:', res.x[0])
print('T1:', res.x[1])
print('T2:', res.x[2])

#--------------------------------------
# Optimization terminated successfully.
#          Current function value: 9.330137
#          Iterations: 6
#          Function evaluations: 1751
# 
# MRE: 9.330136803293637
# k: 1.0239929423715524
# T1: 20.242199361067595
# T2: 10.38718689811098

x_t = [next(pg) for i in mtl]
t, y_A2 = ct.forced_response(A2, T=mtl, U=x_t)

ln = max(abs(max(y_A2)), abs(min(y_A2))) * 1.0

y_A2_g = [y_A2[i]+gauss()*ln for i in range(len(mtl))]

j = 100

max_v = max(y_A2_g[j-4],
            y_A2_g[j-3],
            y_A2_g[j-2],
            y_A2_g[j-1])

min_v = min(y_A2_g[j-4],
            y_A2_g[j-3],
            y_A2_g[j-2],
            y_A2_g[j-1])

mean_v = (min_v + max_v) / 2

res = minimize(f_A2_g,
               [0.01, 0.01, 0.01, mean_v, mean_v, mean_v, mean_v],
               args=(x_t[j:], y_A2_g[j:]),
               method="Powell",
               options={'xtol': 1e-12, 'disp': True},
               bounds=[(0.01, 5), (0.01, 30), (0.01, 20), (min_v, max_v), (min_v, max_v), (min_v, max_v), (min_v, max_v)])

print('\nMRE:', res.fun)
print('k:', res.x[0])
print('T1:', res.x[1])
print('T2:', res.x[2])

#--------------------------------------
# Optimization terminated successfully.
#          Current function value: 10.792331
#          Iterations: 8
#          Function evaluations: 2412
# 
# MRE: 10.792331083024271
# k: 0.9866685845822787
# T1: 19.643955217189404
# T2: 7.761183216604419

Полученные результаты вычисления значений параметров А2-звена вполне удовлетворительны даже при сильном зашумлении данных.

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

Итоги

Основной посыл статьи -- демонстрация возможности детерминирования динамических свойств технологических параметров, описываемых решётчатыми функциями, и возможности прогнозирования изменения значений этих параметров, применяя простые, а оттого высокопроизводительные вычислительные инструменты.

Эти возможности позволят решить, например, следующие практические задачи:

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

Использованная литература

Ротач, В. Я. Теория автоматического управления теплоэнергетическими процессами: Учебник для вузов. — М.: Энергоатомиздат. 1985. — 296 с., ил.
Киреев В. И., Пантелеев А. В. Численные методы в примерах и задачах: Учебное пособие. — 4-е изд., испр. — СПб.: Издательство «Лань», 2015. — 448 с.: ил. — (Учебники для вузов. Специальная литература).
Калиткин Н. Н., Численные методы: учеб. пособие. — 2-е изд., исправленное. — СПб.: БХВ-Петербург, 2011. — 592 с.: ил. — (Учебная литература для вузов)
Березин, И. С., Жидков Н. П. Методы вычислений. — 2-е изд. — М.: Физматлит, 1962. — Т. I.

Линейная регрессия, решётчатые функции и типовые динамические звенья