Зачем это нужно
Метод гармонической линеаризации широко используется для анализа нелинейных систем [1]. Этот метод используется для определения условий возникновения автоколебаний в системах второго и более высокого порядка. При гармонической линеаризации должны выполняться два следующих условия. Замкнутая линейная система должна состоять из двух частей ─ линейной и нелинейной. Линейная часть должна обладать хорошими фильтрующими свойствами для высших гармоник[2]. Системы автоматического управления и регулирования содержат исполнительные механизмы, содержащие нелинейные элементы, поэтому их анализ является весьма актуальной проблемой.
Основной алгоритм
Пусть на вход нелинейного элемента поступает гармонический сигнал
, где 
Преобразовав выходной сигнал в ряд Фурье только для первой гармоники. Получим: 
С учётом (1),(2) общее уравнение нелинейного элемента можно представить в виде[3]:
где ─ действительная «а» – мнимая составляющая коэффициента гармонической линеаризации полученных из условия последующей фильтрации высших гармоник; ─ круговая частота; — оператор дифференцирования.
Постановка задачи
Будем сравнивать два нелинейных элемента приведённых на рис.1,2, с целью выяснить стабильность амплитуды.
Рис.1 Нелинейность типа «насыщение»
Используя (1), пока только для действительной части, получим:
Рис.2 Нелинейность типа «обрыв»
Аналогично используя (1), пока только для действительной части, получим:
Реализация (4),(5) на Python
#!/usr/bin/python
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt # модуль для построения графика
import matplotlib as mpl #модуль для установки шрифта
mpl.rcParams['font.family'] = 'fantasy'# шрифт
mpl.rcParams['font.fantasy'] = 'Comic Sans MS, Aria'# шрифт
import scipy.integrate as spint # модуля для числового интегрирования
import math # модуль для элементарных функций
import numpy as np# модуль для работы с массивами
q=[]# список для коэффициента -насыщение
p=[]#список для коэффициента -обрыв
c=[]# список для границы линейной области
k=1#переменная для установки линейной области
for m in np.arange(1.0,4,0.1): # задание диапазона относительной амплитуды
def f(x): # вычисление интегралов для коэффициентов ряда Фурье
return k*(math.sin(x))**2
a = 0# нижний предел интегрирования
b =math.asin(1/m)# верхний предел интегрирования
z=spint.quad(f,a,b)[0]# результат интегрирования
a =math.pi-math.asin(1/m)#нижний предел интегрирования
b =math.pi#верхний предел интегрирования
z1=spint.quad(f,a,b)[0]#результат интегрирования
q.append(round(z+z1,3))# числовые значения коэффициента - обрыв
def f1(x):# вычисление интегралов для коэффициентов ряда Фурье
return (k/m)*(math.sin(x))
a =math.asin(1/m)#нижний предел интегрирования
b =math.pi-math.asin(1/m)#верхний предел интегрирования
z2=spint.quad(f1,a,b)[0] #результат интегрирования
p.append(round(z+z1+z2,3))#числовые значения коэффициента - насыщение
x=np.arange(1.0,4,0.1)# числовые значения для координаты х для коэффициентов
x2=np.arange(0.0,1.1,0.1)#числовые значения для координаты х для линейной области
for i in np.arange(0.0,1.1,0.1):# диапазон значений для линейной области
c.append(p[0])
plt.title('Коэффициенты гармонической линеаризации', size=14)# заголовок графика
plt.xlabel('Относительная амплитуда сигнала', size=14)# подпись под осью х
plt.ylabel('Приведенный коэффициен усиления', size=14)# подпись под осью у
plt.plot(x, q, color='r', linewidth=2, label='Нелинейность-- обрыв')# построение графика и легенда
plt.plot(x, p, color='b',linewidth=2, label='Нелинейность насыщение')#построение графика и легенда
plt.plot(x2, c, color='g', linewidth=2, label='Линейная область')#построение графика и легенда
plt.legend(loc='best')# лучшее расположение легенды
plt.grid(True)# сетка
plt.show()# визуализация
Результат (да в jpeg, но это график!)
Кстати интегрирование при помощи import scipy.integrate as spint осуществляться стабильно и без ошибок.
Из графика с полной очевидностью следует, что нелинейный элемент, приведенный на рис.2, при условии фильтрации высших гармоник обеспечит большую стабильность амплитуды, чем приведённый на рис.1. Более детальный анализ покажет и меньший процент гармоник, но это тема другой статьи.
Вы, конечно, спросите «А как же другие нелинейности, может они ещё лучше?». Я делал это для низкочастотного (до 5кгц) усилителя специального назначения с типовыми характеристиками, приведенными на рис.1,2. Кстати, характеристики рис.2 нет ни в одном справочнике.
Вариант сокращения кода программы
Можно использовать символьные вычисления, например, так:
from sympy import *
import matplotlib.pyplot as plt
import numpy as np
x=Symbol('x')
m=Symbol('m')
z=integrate(sin(x)**2, (x,0,asinh(m)))+integrate(sin(x)**2, (x,pi-asinh(m),pi)) #символьное преобразование для коэффициент - обрыв
z1=integrate(sin(x)**2, (x,0,asinh(m)))+integrate(m*sin(x), (x,asinh(m),pi-asinh(m)) )+integrate(sin(x)**2, (x,pi-asinh(m),pi))# то же для насыщения
print(z)
print(z1)
Получим:-sin(asinh(m))*cos(asinh(m)) + asinh(m)
2*m*cos(asinh(m)) — sin(asinh(m))*cos(asinh(m)) + asinh(m)
Программное упрощение (simplify(z)) ничего не даёт, поэтому используем соотношения:
, получим: 
Вводим отношение амплитуды к порогу ограничения и окончательно получим:
Соотношение (6) решение уравнения (5), система (7) решение уравнения (4).
Выводы
В статье рассмотрена методика сравнительного анализа нелинейных элементов методом гармонической линеаризации.
Ссылки
- Воронов А.А. Основы теории автоматического управления. Особые
линейные и нелинейные системы. – М.: Энергоиздат, 1981, 300 с. - Гольдфарб Л.С., Балтрушевич А.В., Нетушил А.В. Теория
автоматического управления. – М.: Высшая школа, 1976, 430 с - Хлыпало Е.Н. Нелинейные системы автоматического регулирования / Е.Н. Хлыпало. – Л.: Энергия, 1967. – 343 с.
