31.07.2024 11:16
56h 0 618 Источник

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

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

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

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

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

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

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

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

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

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

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

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

\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 и может быть смело использована для оценки качества прогноза.

Hidden text
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_2T_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<1T^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

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

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_i3\,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_iy_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_iy_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

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

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_iy_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:]

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

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

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

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

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_i11\,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_iy_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 рассчитываются как комбинации параметров звена, в случае A2-звена — 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.

Hidden text

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})
Hidden text
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, которые для поиска глобального минимума функции многих переменных используют нейросетевые инструменты.

Итоги

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

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

  • реализация алгоритма автоматической настройки регуляторов систем автоматического регулирования по известным динамическим свойствам объекта управления;

  • реализация алгоритма предиктивного ПИД-регулятора, в котором наряду с П-, И-, Д- компонентами присутствует компонента прогноза изменения регулируемого технологического параметра;

  • создание модели технологического многопараметрического объекта управления, работающей в режиме реального времени, которая позволяет решить следующие задачи:

    • проведение исследований и экспериментов над технологическим объектом управления при разработке новых решений в управлении;

    • создание полнофункциональных тренажёрных комплексов систем управления реальными технологическими объектами;

    • создание систем прогнозирования хода технологического процесса и формирование рекомендаций оператору по безопасному и более эффективному его ведению.

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

  1. Ротач, В. Я. Теория автоматического управления теплоэнергетическими процессами: Учебник для вузов. — М.: Энергоатомиздат. 1985. — 296 с., ил.

  2. Киреев В. И., Пантелеев А. В. Численные методы в примерах и задачах: Учебное пособие. — 4-е изд., испр. — СПб.: Издательство «Лань», 2015. — 448 с.: ил. — (Учебники для вузов. Специальная литература).

  3. Калиткин Н. Н., Численные методы: учеб. пособие. — 2-е изд., исправленное. — СПб.: БХВ-Петербург, 2011. — 592 с.: ил. — (Учебная литература для вузов)

  4. Березин, И. С., Жидков Н. П. Методы вычислений. — 2-е изд. — М.: Физматлит, 1962. — Т. I.

Комментарии (0)