Оценка параметров модели электропривода

Часто при моделировании возникает ситуация, когда законы, описывающие поведение объекта, известны, однако некоторые параметры, входящие в уравнения, остаются неопределёнными. В данном примере поэтапно рассмотрен процесс оценки неизвестных параметров модели электропривода на основе экспериментальных данных.

Подключение необходимых пакетов

Установим необходимые пакеты:

import Pkg
Pkg.add(["Interpolations","NLopt"])

Подключим необходимые пакеты:

using DataFrames, Plots,Interpolations, NLopt, MATLAB

Модель электропривода

Рассмотрим модель электропривода "motor_pe", схема которой представлена на рисунке ниже. Она состоит из контроллера, усилителя и электродвигателя с нагрузкой. Вход модели представляет собой сигнал управляющий положением нагрузки, выход - положение нагрузки.

Для запуска симуляции модели нам необходимо указать значения следующих параметров:

J1 - инерция ротора электродвигателя, ;
J2 - инерция нагрузки, ;
k - жесткость гибкого вала, ;
Km - постоянная электродвигателя, ;
Ka - коэффициент усиления усилителя, ;
b1 - коэффициент демпфирования подшипников электродвигателя, ;
b2 - коэффициент демпфирования подшипников нагрузки, ;
b12 - коэффициент демпфирования гибкого вала, .

Некоторые из этих параметров известны и могут быть получены из технической документации или специализированных программ (например, систем автоматизированного проектирования — САПР), включая коэффициент усиления усилителя, моменты инерции нагрузки и ротора электродвигателя. Другие параметры могут быть определены с помощью прямых измерений, например, жесткость гибкого вала. Введём переменные, содержащие значения известных параметров:

J1 = 1.00e-6; 
J2 = 1.35e-7;    
k = 0.001;      
Km = 16.4e-3;   
Ka = 1.3/2.5;

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

b1 = 1e-5;      
b2 = 1e-5;     
b12 = 1e-5;

Загрузка результатов измерений

Загрузим значения измеренных сигналов из файла с экспериментальными данными "motor_data". В вектор tdata запишем значения моментов времени, в которые проводились измерения, в вектор udata — входной сигнал, представляющий заданное положение нагрузки, а в вектор ydata — выходной сигнал, соответствующий измеренному положению нагрузки:

data=read_matfile("$(@__DIR__)/motor_data.mat")
pa_input_voltage = jvector(data["pa_input_voltage"])
pa_output_voltage = jvector(data["pa_output_voltage"])
tdata = jvector(data["tdata"]) 
udata = jvector(data["udata"]) 
ydata = jvector(data["ydata"]);

Построим графики измеренных сигналов:

gr()
input = plot(tdata, udata, title = "Вход", xlabel="t, c",legend=false)
output_measure = plot(tdata, ydata, title = "Выход", xlabel="t, c", legend=false)
plot(output_measure,input,layout = (2, 1))

Сравнение результатов моделирования и измерений

Для сравнения результатов моделирования и экспериментальных измерений необходимо подать входной сигнал в модель и зарегистрировать её выходные данные.
Создадим переменную типа WorkspaceArray для передачи входного сигнала udata в модель. Следует отметить, что в качестве значений времени для входного сигнала необходимо использовать вектор tdata.

udata_workspace = WorkspaceArray("udata_workspace", DataFrame(time = tdata, value = udata));

С помощью блока Заданное положение нагрузки (выделенного жёлтым на изображении ниже) входной сигнал udata_workspace передаётся на вход контроллера. Выходной сигнал load_position регистрируется при помощи логирования, его значения становятся доступны после завершения симуляции. Блок Заданное положение нагрузки и сигнал load_position соответствуют точкам сбора входных и выходных данных реальной системы.

Откроем модель "motor_pe", выполним её симуляцию, а затем построим график, отображающий результаты эксперимента и моделирования:

model = engee.load("motor_pe.engee")
logout = engee.run(model, verbose=false)
output = plot(tdata, ydata, title = "Выход", xlabel="t, c", label = "эксперимент", legend=true)
plot!(logout["load_position"].time, logout["load_position"].value, label = "моделирование", linewidth=2)
plot(output,input,layout = (2, 1))

Как видно из графика, динамика реальной системы и модели не совпадают. Это указывает на необходимость более точной настройки неизвестных параметров b1, b2, b12.

Оценка параметров при помощи оптимизации

Необходимо определить такие значения параметров b1, b2, b21, при которых выходной сигнал модели, записываемый блоком Положение нагрузки, максимально точно соответствует измеренным данным ydata.

Для этого будем минимизировать сумму квадратов ошибок между измеренным и смоделированным выходными сигналами. Ошибки вычисляются как разность между значениями измеренного выходного сигнала ydata и выходного сигнала модели load_position в моменты времени, заданные вектором tdata. Для минимизации суммы квадратов ошибок используем один из доступных алгоритмов оптимизации.

Важно отметить, что оцениваемые параметры b1, b2, b21 должны быть переменными (отображаться в рабочей области) и использоваться в модели. В данном примере они задействованы в блоке motor_pe/Электродвигатель и нагрузка/State-Space.

# @markdown # Оценка параметров модели
# @markdown ## Для начала оценки параметров модели введите:
# @markdown ### Название модели:
# Открываем модель
model_name = "motor_pe" # @param {type:"string"}
if model_name in [m.name for m in engee.get_all_models()] # Проверка условия загрузки модели в ядро
    model = engee.open(model_name) # Открыть модель
    model = engee.gcm() 
else
    @error "Модель $model_name не открыта."
end
# @markdown ### Вход (блок типа FromWorkspace):
input_name = "Заданное положение нагрузки" # @param {type:"string"}
input_path = engee.find_system(model; depth=0, blockparams=["BlockType" => "FromWorkspace","BlockName"=>input_name])
# @markdown ### Значение входного сигнала:
input_data = [tdata udata] # @param {type:"raw"}
typeof(input_data) === Matrix{Float64} || throw(ArgumentError("Входной сигнал должен быть матрицей содержащей числа типа Float64."))
size(input_data,2) == 2 || throw(ArgumentError("Входной сигнал должен иметь 2 столбца - время и значения сигнала."))

isempty(engee.find_system(model; depth=0, blockparams=["BlockType" => "FromWorkspace","BlockName"=>input_name])) && throw(ArgumentError("Блока $input_name не существует в корневой системе модели $(model.name)."))
pe_udata_workspace = WorkspaceArray("pe_udata_workspace", DataFrame(time = input_data[:,1], value = input_data[:,2]));
engee.set_param!("$(model.name)/$input_name", "VariableName" =>"pe_udata_workspace") 

# @markdown ### Выход (логируемый сигнал):
output_name = "load_position" # @param {type:"string"}
output_name in keys(engee.run(model, verbose=false)) || throw(ArgumentError("Сигнал $output_name не найден среди логируемых сигналов."))

# @markdown ### Значение выходного сигнала:
output_data = [tdata ydata] # @param {type:"raw"}
typeof(output_data) === Matrix{Float64} || throw(ArgumentError("Выходной сигнал должен быть матрицей, содержащей числа типа Float64."))
size(output_data,2) == 2 || throw(ArgumentError("Выходной сигнал должен иметь 2 колонки - время и значения сигнала."))

# @markdown ### Имена оцениваемых параметров:
# Создаем вектор переменных, их начальные значения и ограничения
p_names = ["b1","b2","b21"] # @param {type:"raw"}
p_names isa Vector{String} || throw(ArgumentError("Названия параметров должны быть введены в виде вектора строк."))
parametrs_number = length(p_names)
parametrs_number>0 || throw(ArgumentError("Количество параметров должно быть больше нуля."))

# @markdown ### Начальные значения параметров:
p_0 = [1e-5, 1e-5, 1e-5] # @param {type:"raw"}
p_0 isa Vector{Float64} || throw(ArgumentError("Начальные значения парметров должны быть вектором, содержащим числа типа Float64."))
parametrs_number==length(p_0) || throw(ArgumentError("Количество начальных значений и парметров должно совпадать."))

# @markdown ### Нижняя граница значений параметров:
p_lower_bounds = [0.0, 0.0, 0.0] # @param {type:"raw"}
p_lower_bounds isa Vector{Float64} || throw(ArgumentError("Нижняя граница значений параметров должна быть вектором, содержащим числа типа Float64."))
parametrs_number==length(p_lower_bounds) || throw(ArgumentError("Количество значений нижней границы и парметров должно совпадать."))

# @markdown ### Верхняя граница значений параметров:
p_upper_bounds = [1e-3,1e-3,1e-3] # @param {type:"raw"}
p_upper_bounds isa Vector{Float64} || throw(ArgumentError("Верхняя граница значений параметров должна быть вектором, содержащим числа типа Float64."))
parametrs_number==length(p_0) || throw(ArgumentError("Количество начальных значений и парметров должно совпадать."))

# @markdown ## Настройки оптимизации:
# @markdown ### Алгоритм оптимизации ("G" в начале названия алгоритма говорит о том что он глобальный, "L" - локальный)
algoritm = :G_MLSL # @param [:GN_DIRECT,:GN_DIRECT_L,:GN_DIRECT_L_RAND,:GN_CRS2_LM,:G_MLSL,:GN_AGS,:GN_ISRES,:GN_ESCH,:LN_COBYLA,:LN_BOBYQA,:LN_PRAXIS,:LN_NELDERMEAD,:LN_SBPLX] {type:"raw"}
# @markdown ### Вспомогательный алгоритм (только для :G_MLSL)
local_algoritm = :LN_COBYLA # @param [:LN_COBYLA,:LN_BOBYQA,:LN_PRAXIS,:LN_NELDERMEAD,:LN_SBPLX] {type:"raw"}
# Настройка параметров оптимизации
# @markdown ### Минимальное относительное изменение шага
xtol_rel = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("Минимальное относительное изменение шага должно быть положительным."))
# @markdown ### Минимальное абсолютное изменение шага
xtol_abs = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("Минимальное абсолютное изменение шага должно быть больше положительным."))
# @markdown ### Минимальное относительное изменение целевой функции
ftol_rel = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("Минимальное относительное изменение целевой функции должно быть положительным."))
# @markdown ### Минимальное абсолютное изменение целевой функции
ftol_abs = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("Минимальное абсолютное изменение целевой функции должно быть положительным."))
# @markdown ### Максимальное количество оценок целевой функции
maxeval = 30 # @param {type:"integer"}
xtol_rel>0 || throw(ArgumentError("Максимальное количество оценок целевой функции быть больше нуля."))

# @markdown ## Отображение результатов оценки параметров
# @markdown ### График измеренных данных и результатов моделирования:
compare = true # @param {type:"boolean"}
# @markdown ### График ошибок между измеренными данными и результатами моделирования:
reside = true # @param {type:"boolean"}
# @markdown ### График значений оценок параметров:
variables = true # @param {type:"boolean"}
# @markdown ### График значений целевой функции:
objective = true # @param {type:"boolean"}

#График сравнения выхода модели и измеренных данных
function compare_plot(model, output_name::String, tdata::Vector{Float64}, udata::Vector{Float64}, ydata::Vector{Float64})
    logout = engee.run(model, verbose=false)
    out_plot = plot(
        logout[output_name].time, 
        logout[output_name].value, 
        label="моделирование", 
        title = "Выход", 
        xlabel="t, с", 
        legend=:bottomright,
        linewidth=3
        )
    plot!(out_plot, tdata, ydata, label="эксперимент")
    in_plot = plot(tdata, udata, title = "Вход", xlabel="t, c", legend = false)
    p = plot(out_plot,in_plot,layout = (2, 1), plot_title="Эксперимент и моделирование")
end

# график невязок
function residual_plot(model, output_name::String, tdata::Vector{Float64}, ydata::Vector{Float64})
    logout = engee.run(model, verbose=false)
    nearest = interpolate((logout[output_name].time,), logout[output_name].value, Gridded(Constant()))
    residals = ydata .- nearest.(tdata)
    reside = plot(tdata, residals, title = "Ошибки", xlabel="t, с", legend = false)
end

# График измения переменых
function variable_plot(p_names::Vector{String}, parametrs_values)
    p_trace = plot(title = "Значения параметров", xlabel="Количество оценок целевой функции",)
    parametrs_values_m = mapreduce(permutedims, vcat, parametrs_values)
    for i in 1:n
        plot!(p_trace, 1:length(parametrs_values_m[:,i]), parametrs_values_m[:,i], markers=:auto, label = p_names[i])
    end
    return p_trace
end

n = parametrs_number

# Векторы для логирования значений целевой функции и переменных
loss_function_vector = Vector{Float64}()
parametrs_values = Vector{Vector{Float64}}()

#функция потерь
function myfunc(x::Vector, grad::Vector, output_name::String, measure_time, measure_value, parametrs_name...)
    str = ""
    iterations_value = Vector{Float64}()
    for (value,p) in zip(x, parametrs_name)
        global parametr = Symbol(p)
        @eval(($parametr)=$(value))
        str = str*"$p = $value, "
        push!(iterations_value, value)
    end
    push!(parametrs_values, iterations_value)
    simulation_out = engee.run(model, verbose=false) |> x->collect(x[output_name])
    if measure_time[end]>simulation_out.time[end]  
        push!(simulation_out,[measure_time[end], simulation_out.time[end]])
    end
    nearest = interpolate(
        (simulation_out.time,), 
        simulation_out.value, 
        Gridded(Constant())
    )
    l = sum(abs2, measure_value .- nearest.(measure_time))
    str = "$(length(parametrs_values)). "*str*"loss = $l"
    println(str)
    push!(loss_function_vector, l)
    return l
end


opt = Opt(algoritm, n) # алгоритм без градиента
algoritm === :G_MLSL && NLopt.local_optimizer!(opt, Opt(local_algoritm, n))
opt.lower_bounds = p_lower_bounds
opt.upper_bounds = p_upper_bounds
opt.xtol_rel = xtol_rel
opt.xtol_abs = xtol_abs
opt.ftol_rel = ftol_rel
opt.ftol_abs = ftol_rel
opt.maxeval =  maxeval

opt.min_objective = (x,g) -> myfunc(x,g,output_name,output_data[:,1],output_data[:,2],p_names...)

(minf,minx,ret) = optimize(opt,p_0)
numevals = opt.numevals 

str = "\nРезультаты оценки:\n"
for (name, value) in zip(p_names, minx) 
    str = str*"$name = $value\n"
end
str = str*"\nЗначение целевой функции = $minf.\n\nОптимизация завершена после $numevals оценок целевой функции cо статусом: $ret.\n"
print(str)

compare && display(compare_plot(model, output_name, input_data[:,1], input_data[:,2], output_data[:,2]))
println("")
reside && display(residual_plot(model, output_name, output_data[:,1], output_data[:,2]))
println("")
objective &&
    display(
        plot(1:numevals, loss_function_vector, title = "Целевая функция", xlabel="Количество оценок целевой функции", legend = false, markers = :circle)
    )
println("")   
variables && display(variable_plot(p_names, parametrs_values))
println("")

1. b1 = 1.0e-5, b2 = 1.0e-5, b21 = 1.0e-5, loss = 131163.2887540543
2. b1 = 0.00042129156796609637, b2 = 0.0001388549111491948, b21 = 0.00026820318603994474, loss = 603468.7780001939
3. b1 = 0.0003553762405112313, b2 = 0.0003798108457121816, b21 = 0.000512887209368987, loss = 614893.4957088778
4. b1 = 0.0008082835815970666, b2 = 5.6455325296927654e-5, b21 = 0.0006708533954141649, loss = 605100.485177848
5. b1 = 0.00019251590525647854, b2 = 0.0009458991230097862, b21 = 0.0006315424151859726, loss = 615394.4066690244
6. b1 = 1.0e-5, b2 = 1.0e-5, b21 = 1.0e-5, loss = 131163.2887540543
7. b1 = 1.7500000000000002e-5, b2 = 1.0e-5, b21 = 1.0e-5, loss = 209221.61878551292
8. b1 = 1.0e-5, b2 = 1.7500000000000002e-5, b21 = 1.0e-5, loss = 208161.7876511638
9. b1 = 1.0e-5, b2 = 1.0e-5, b21 = 1.7500000000000002e-5, loss = 131163.2887540543
10. b1 = 4.660575198320843e-6, b2 = 4.733070838978989e-6, b21 = 1.0e-5, loss = 34003.99464283905
11. b1 = 8.470329472543003e-22, b2 = 1.6940658945086007e-21, b21 = 1.0e-5, loss = 26592.070070605747
12. b1 = 3.2487355184912022e-6, b2 = 2.0816681711721759e-22, b21 = 8.126976366701769e-6, loss = 5318.483821111458
13. b1 = 3.159138766157072e-6, b2 = 1.6132928326584313e-21, b21 = 4.3780468632339875e-6, loss = 5213.718258001723
14. b1 = 6.907644735263423e-6, b2 = 1.4571677198205197e-21, b21 = 4.272202640084758e-6, loss = 18773.09903792265
15. b1 = 1.2843838515077049e-6, b2 = 1.128474576789396e-36, b21 = 4.408362040918534e-6, loss = 11211.56421356133
16. b1 = 3.1591387661570723e-6, b2 = 6.314226320581347e-7, b21 = 4.3780468632339875e-6, loss = 5854.704219020841
17. b1 = 5.0331408497233005e-6, b2 = 0.0, b21 = 4.3168815888584465e-6, loss = 9721.443209286368
18. b1 = 4.096139807940186e-6, b2 = 0.0, b21 = 4.347464226046217e-6, loss = 6605.131199555274
19. b1 = 3.1558740030327797e-6, b2 = 1.6143061990385022e-21, b21 = 4.277576733299245e-6, loss = 5283.599069547135
20. b1 = 2.323527217024245e-6, b2 = 0.0, b21 = 4.803087554175263e-6, loss = 6096.148521317298
21. b1 = 2.880871606471125e-6, b2 = 2.6020852139651799e-23, b21 = 3.4827962626531693e-6, loss = 5241.365644792966
22. b1 = 3.6108253578913063e-6, b2 = 0.0, b21 = 4.25272413844909e-6, loss = 5663.868162772741
23. b1 = 2.937563588077814e-6, b2 = 3.2265856653168625e-21, b21 = 4.4544411109503256e-6, loss = 5274.899960128769
24. b1 = 3.1860433362031623e-6, b2 = 1.6116526510050095e-21, b21 = 4.406208499303296e-6, loss = 5247.156054409827
25. b1 = 3.137707490175435e-6, b2 = 2.6020852139651987e-23, b21 = 4.144653759148201e-6, loss = 5256.7638593533975
26. b1 = 3.1591387661570714e-6, b2 = 1.0336026984674575e-7, b21 = 4.3780468632339875e-6, loss = 5305.595047753606
27. b1 = 3.371394146001847e-6, b2 = 3.2265856653168625e-21, b21 = 4.477441502686777e-6, loss = 5375.470859753354
28. b1 = 3.2652664560794595e-6, b2 = 3.2265856653168625e-21, b21 = 4.427744182960382e-6, loss = 5283.371686285544
29. b1 = 3.103487872268651e-6, b2 = 0.0, b21 = 4.396382771951099e-6, loss = 5253.502820320646
30. b1 = 3.1641847304454216e-6, b2 = 3.2265856653168625e-21, b21 = 4.319670791222732e-6, loss = 5267.258290082376

Результаты оценки:
b1 = 3.159138766157072e-6
b2 = 1.6132928326584313e-21
b21 = 4.3780468632339875e-6

Значение целевой функции = 5213.718258001723.

Оптимизация завершена после 30 оценок целевой функции cо статусом: MAXEVAL_REACHED.

После выполнения необходимого количества итераций были получены следующие значения оцениваемых параметров:

b1 = 2.9731e-6 — коэффициент демпфирования подшипников электродвигателя;
b2 = 9.1799e-23 — коэффициент демпфирования подшипников нагрузки;
b12 = 7.439e-6 — коэффициент демпфирования гибкого вала.

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

Качество оценки параметров зависит от следующих факторов:

точности начальных значений параметров;
наложенных ограничений на значения параметров;
настроек алгоритма оптимизации;
количества оцениваемых параметров.

Результаты

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

Оценка параметров модели электропривода