Проектирование САУ фокусировки камеры смартфона

Подключим необходимые библиотеки.

import Pkg; 
Pkg.add(["ControlSystemIdentification", "ControlSystems"])

using ControlSystems                # основная библиотека для работы с САУ
using ControlSystemIdentification   # библиотека для оценки динамических моделей на основе данных 

Объект управления

Модель, приведенная ниже, симулирует систему управления автофокусом камеры смартфона.

Рассмотрим модель автофокусировки смартфона. Электроника смартфона осуществляет управление мотором за счет обратной связи от камеры. Абстрагируемся от технологии и скажем, что нам нужно чтобы переходный процесс занимал не более 0.6 секунд с момента получения команды на фокусировку. С помощью ПИД регулятора попробуем реализовать такую систему.

Наша система мотор-линза является нелинейной, поэтому прежде всего линеаризуем модель.

Линеаризация объекта управления

Запустим модель с помощью командного управленния. В результате в рабочей области получим переменную data, которая содержит результаты симуляции.

modelName = "Линза_ОУ"
if !(modelName in [m.name for m in engee.get_all_models()]) engee.load( "$(@__DIR__)/$(modelName).engee"); end;
data = engee.run( modelName )

SimulationResult(
    "Преобразование единиц измерения.1" => WorkspaceArray{Float64}("Линза_ОУ/Преобразование единиц измерения.1")
,
    "Ступенчатая функция.1" => WorkspaceArray{Float64}("Линза_ОУ/Ступенчатая функция.1")

)

Запишем в переменные данные входных и выходных сигналов в переменные и построим их на графике.

dataIn = collect(data["Ступенчатая функция.1"]).value;
dataOut = collect(data["Преобразование единиц измерения.1"]).value;
data_t = collect(data["Преобразование единиц измерения.1"]).time;

plot( data_t, [dataIn dataOut], label=["Входной сигнал" "Выходной сигнал"] )

Перейдем к линеаризации модели. Все методы для идентификации модели библиотеки ControlSystemIdentification принимают на вход объект типа AbstractIdData. Получить объект идентификационных данных временной области позволяет функция iddata. Далее используя метод subspaceid получим описание системы в виде пространства состояний.

Ts = data_t[2] - data_t[1]                     # Получим частоту дискретизации системы из вектора времени
d  = iddata(dataOut, dataIn, Ts);              # создаем объект идентификации данных во временном пространстве 
sys_d = subspaceid(d, :auto);                  # запускаем метод идентификации модели пространства состояний

Сравним полученное описание модели и с моделью собранной в среде моделирования. Для этого используем блок Пространство состояний и в качестве значений матриц запишем туда полученные значения A,B,C,D.

sys = d2c(sys_d);   # перейдем к непрерывной модели
A = sys.A;          # находим матрицу A
B = sys.B;          # находим матрицу B
C = sys.C;          # находим матрицу C
D = sys.D;          # находим матрицу D

┌ Warning: d2c does not handle a non-zero cross covariance S, S will be set to zero
└ @ ControlSystemIdentification /usr/local/ijulia-demos/packages/ControlSystemIdentification/KM4Nc/src/ControlSystemIdentification.jl:359

Подбор коэффициентов: функция pid()

Чтобы оценить насколько подходят подобранные значения коэффициентов, будем строить графики ЛАЧХ, ЛФЧХ и годограф Найквиста. А также оценивать качество переходного процесса. Для удобного отображения графиков частотных характеристик ниже реализована функция pid_marginplot. Она подставляет значения полученной системы у правления в функции для построения частотных характеристик: marginplot и nyquistplot, и выводит их на графики.

function pid_marginplot(C)
    f1 = marginplot(sys_tf * C, w)
    f2 = nyquistplot(sys_tf * C, xlims=(-10,1), ylims=(-2, 10), title="")
    plot(f1,f2, titlefontsize=10)
end

pid_marginplot (generic function with 1 method)

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

sys_tf = tf(sys);
w = exp10.(LinRange(-2.5, 3, 500));

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

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

kp = 0.4 # @param {type:"slider", min:0.1, max:5, step:0.1}
ki = 1.322 # @param {type:"slider", min:0.0, max:20, step:0.001}
kd = 0.03 # @param {type:"slider", min:0.0, max:0.5, step:0.01}
Tf = 0.012 # @param {type:"slider", min:0.0, max:1, step:0.001}

C0 = pid(kp, ki, kd; Tf, form=:parallel, filter_order=2)
pid_marginplot(C0)

Проверим устойчива ли наша система. Это позволяет сделать функция isstable.

isstable(feedback(sys_tf * C0))     # Проверка: устойчива ли система с полученным ПИД регулятором

true

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

plot(stepinfo(step(feedback(sys_tf * C0), 5)))      # Оценка переходного процесса САУ

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

Подбор коэффициентов: функция loopshapingPID()

Подбор коэффициентов с помощью функции loopshapingPID основан на частотном подходе (Loop Shaping), который позволяет задать желаемую форму АЧХ разомкнутой системы для обеспечения нужных характеристик замкнутого контура. Таким образом это метод удобен, если известны желаемые частотные свойства системы (например, требуемая полоса пропускания или уровень подавления возмущений).

Попробуем применить данный метод к нашему примеру и получить желаемые характеристики.

Функция принимает описание системы в виде передаточной функции и частоту, на которой желательно достичь единичное усиление разомкнутой системы. Дополнительно функция принимает значения параметров дополнительной функции чувствительности Mt и желаемый фазовый сдвиг ϕt. По умолчанию эти параметры имеют оптимальные значения, но их можно изменить. Также можно настроить форму синтезируемого ПИД регулятора и вывести графики, которые помогут оценить качество полученной САУ.

ω  = 17
Tf = 1/1000ω   
C0, kp, ki, kd, fig, CF = loopshapingPID(sys_tf, ω; Mt = 1.3, ϕt = 75, form=:parallel, doplot = true)

(C = TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
0.0493284587237057s^2 + 1.0063456913621145s + 6.26737143322893
--------------------------------------------------------------
                             1.0s

Continuous-time transfer function model, kp = 1.0063456913621145, ki = 6.26737143322893, kd = 0.0493284587237057, fig = Plot{Plots.PlotlyJSBackend() n=18}, CF = TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
0.0493284587237057s^2 + 1.0063456913621145s + 6.26737143322893
--------------------------------------------------------------
   3.4602076124567474e-9s^3 + 8.31890330807703e-5s^2 + 1.0s

Continuous-time transfer function model)

Давайте попробуем изменить частоту. Если анализировать диаграмму Найквиста, то можно заметить, что характеристика проходит очень близко к критической точки. Давайте попробуем скорректировать это и установим ω = 17. Теперь характеристика скорректированной системы будет касаться окружности Mt в точке, где ω = 17 под углом ϕt = 75.

plot(fig)

Убедимся в устойчивости полученной системы.

isstable(feedback(CF*sys_tf))

true

Построим ЛАЧХ и ЛФЧХ системы и отразим запасы устойчивости.

marginplot([sys_tf, CF*sys_tf])

Также постороим переходный процесс полученной САУ. Видим, что время переходного процесса около 0.4 сек, что удовлетворяет нашим требованиям.

plot(stepinfo(step(feedback(CF * sys_tf), 1)))

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

Блоки, использованные в примере

• Математическая функция
• Деление
• Интегратор второго порядка
• Ступенчатая функция
• Сложение