Estimation of missing model parameters

If we think about creating models of a control object in terms of the data we have about the model, then there are often 3 approaches:
A white box model is a model in which we know all the basic dynamics of the system.
On the contrary, the black box model implies that we don't know the dynamics at all.
This example examines the intermediate gray box approach, when we know the structure of the model, but we lack some parameters.
The missing parameters can be selected based on experimental data. This is a rather complex algorithm, but thanks to the collaboration of scripts and models, even such a difficult task can be implemented in the code.

Loading data from a table

using DataFrames, CSV;

Import_data = "/user/start/examples/controls/parameter_estimation_webinar/parameter_estimation.csv" # @param {type:"filepicker"}
table = CSV.read(Import_data, DataFrame);
print(".csv file imported")
Показать_таблицу = true # @param {type:"boolean"}
if Показать_таблицу
    display(table)
end

.csv файл импортирован

Столбец_времени = 1 # @param {type:"integer"}
Столбец_входа = 2 # @param {type:"integer"}
Столбец_выхода = 3 # @param {type:"integer"}
Начинать_со_строки = 1 # @param {type:"integer"}
exp_time_data = table[Начинать_со_строки:end, Столбец_времени]; 
exp_in_data = table[Начинать_со_строки:end, Столбец_входа];
exp_out_data = table[Начинать_со_строки:end, Столбец_выхода];

In this example, the control object is a DC motor. The model contains all the known parameters from the engine documentation. But the friction is not given in the documentation, so now there are approximate values in the blocks of coefficients of viscous friction and sliding friction.

The next cell runs the model for simulation to make sure that the results of the model will be far from the experimental data.

Comparison of experimental data and model data

model_name = "parameter_estimation" # @param {type:"string"}
input_name = "Из рабочей области" # @param {type:"string"}

if model_name in [m.name for m in engee.get_all_models()]
    m = engee.open( model_name ) # loading the model
else
    m = engee.load("$(model_name).engee" )
end

df_in = DataFrame(time=exp_time_data, value=exp_in_data)
df_in = Float64.(df_in);
wa_in = WorkspaceArray(string(rand()), df_in)
engee.set_param!("$(model_name)/$input_name", "VariableName" =>"wa_in") 


logout = engee.run(m, verbose=false)
input = plot(exp_time_data, exp_in_data, title = "Entrance", xlabel="t, c", label = "experiment", legend=true)
output = plot(exp_time_data, exp_out_data, title = "Exit", xlabel="t, c", label = "experiment", legend=true)
plot!(logout["w"].time, logout["w"].value, label = "simulation", linewidth=2)
plot(input, output, layout = (2, 1))

Connecting libraries for evaluation

using DataFrames
try
    using NLopt
catch
    Pkg.add(["NLopt"])
    using NLopt
end
try
    using Interpolations
catch
    Pkg.add(["Interpolations"])
    using Interpolations
end

# @markdown # Estimation of model parameters
# @markdown ## To start evaluating the model parameters, enter:
# @markdown ### Model name:
# Opening the model
model_name = "parameter_estimation" # @param {type:"string"}
if model_name in [m.name for m in engee.get_all_models()] # Checking the condition for loading a model into the kernel
    model = engee.open(model_name) # Open the model
    model = engee.gcm() 
else
    @error "Модель $model_name не открыта."
end
# @markdown ### Login (FromWorkspace type block):
input_name = "Из рабочей области" # @param {type:"string"}
input_path = engee.find_system(model; depth=0, blockparams=["BlockType" => "FromWorkspace","BlockName"=>input_name])
# @markdown ### Input value:
input_data = [exp_time_data exp_in_data] # @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 (logable signal):
output_name = "w" # @param {type:"string"}
output_name in keys(engee.run(model, verbose=false)) || throw(ArgumentError("Сигнал $output_name не найден среди логируемых сигналов."))

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

# @markdown ### Names of evaluated parameters:
# Creating a vector of variables, their initial values and constraints
p_names = ["b","c"] # @param {type:"raw"}
p_names isa Vector{String} || throw(ArgumentError("Названия параметров должны быть введены в виде вектора строк."))
parametrs_number = length(p_names)
parametrs_number>0 || throw(ArgumentError("Количество параметров должно быть больше нуля."))

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

# @markdown ### Lower bound of parameter values:
p_lower_bounds = [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 ### Upper bound of parameter values:
p_upper_bounds = [1.0, 10.0] # @param {type:"raw"}
p_upper_bounds isa Vector{Float64} || throw(ArgumentError("Верхняя граница значений параметров должна быть вектором, содержащим числа типа Float64."))
parametrs_number==length(p_0) || throw(ArgumentError("Количество начальных значений и парметров должно совпадать."))

# @markdown ## Optimization Settings:
# @markdown ### Optimization algorithm (the "G" at the beginning of the algorithm name indicates that it is global, the "L" is local)
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 ### Auxiliary algorithm (only for :G_MLSL)
local_algoritm = :LN_COBYLA # @param [:LN_COBYLA,:LN_BOBYQA,:LN_PRAXIS,:LN_NELDERMEAD,:LN_SBPLX] {type:"raw"}
# Configuring optimization parameters
# @markdown ### Minimum relative step change
xtol_rel = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("Минимальное относительное изменение шага должно быть положительным."))
# @markdown ### Minimum absolute step change
xtol_abs = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("Минимальное абсолютное изменение шага должно быть больше положительным."))
# @markdown ### Minimal relative change in objective function
ftol_rel = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("Минимальное относительное изменение целевой функции должно быть положительным."))
# @markdown ### The minimum absolute change in the objective function
ftol_abs = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("Минимальное абсолютное изменение целевой функции должно быть положительным."))
# @markdown ### Maximum number of estimates of the objective function
maxeval = 30 # @param {type:"integer"}
xtol_rel>0 || throw(ArgumentError("Максимальное количество оценок целевой функции быть больше нуля."))

# @markdown ## Displaying parameter evaluation results
# @markdown ### Graph of measured data and simulation results:
compare = true # @param {type:"boolean"}
# @markdown ### Error graph between measured data and simulation results:
reside = true # @param {type:"boolean"}
# @markdown ### Graph of parameter estimates:
variables = true # @param {type:"boolean"}
# @markdown ### Graph of objective function values:
objective = true # @param {type:"boolean"}

# Graph comparing the output of the model and the measured data
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="simulation", 
        title = "Exit", 
        xlabel="t, with", 
        legend=:bottomright,
        linewidth=3
        )
    plot!(out_plot, tdata, ydata, label="experiment")
    in_plot = plot(tdata, udata, title = "Entrance", xlabel="t, c", legend = false)
    p = plot(out_plot,in_plot,layout = (2, 1), plot_title="Эксперимент и моделирование")
end

# schedule of discrepancies
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 = "Mistakes", xlabel="t, with", legend = false)
end

# Schedule of change of variables
function variable_plot(p_names::Vector{String}, parametrs_values)
    p_trace = plot(title = "Parameter values", xlabel="Number of estimates of the objective function",)
    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

# Vectors for logging the values of the objective function and variables
loss_function_vector = Vector{Float64}()
parametrs_values = Vector{Vector{Float64}}()

# loss function
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) # gradient-free algorithm
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 = "Target function", xlabel="Number of estimates of the objective function", legend = false, markers = :circle)
    )
println("")   
variables && display(variable_plot(p_names, parametrs_values))
println("")

1. b = 0.01, c = 0.7, loss = 782827.2305232552
2. b = 0.6269504073816847, c = 3.165466261913612, loss = 2.0824007980727735e6
3. b = 0.19966846490239398, c = 2.6384859605251734, loss = 1.5577384910763318e6
4. b = 0.43616990062596706, c = 5.260365036287981, loss = 2.1092703100347244e6
5. b = 0.13062205726902043, c = 0.10527521129684247, loss = 491882.1601205184
6. b = 0.13062205726902043, c = 0.10527521129684247, loss = 491882.1601205184
7. b = 0.22858860022078575, c = 0.10527521129684247, loss = 994688.4244879746
8. b = 0.13062205726902043, c = 0.18423161976947433, loss = 511161.069044274
9. b = 0.03272744808240044, c = 0.10225003874684467, loss = 153347.19109891905
10. b = 0.0, c = 0.02782974568791464, loss = 3.8266403721268633e6
11. b = 0.02681806322965815, c = 0.14143990273494417, loss = 290065.37337741273
12. b = 0.0571202535648882, c = 0.10402153357535857, loss = 11104.510447718792
13. b = 0.05772050887137111, c = 0.10013441411262017, loss = 12468.798951365312
14. b = 0.08125046811420773, c = 0.10740009755882907, loss = 140047.86862450984
15. b = 0.05638690894821858, c = 0.1081700517083551, loss = 9544.510125887842
16. b = 0.08042415294727309, c = 0.11195469496423044, loss = 135148.34762457936
17. b = 0.05502463241104411, c = 0.11520759135488587, loss = 7055.968949365142
18. b = 0.030655512945794418, c = 0.11323568367800015, loss = 194849.15061134959
19. b = 0.042840072678419265, c = 0.114221637516443, loss = 27108.42978479694
20. b = 0.06663619406442059, c = 0.11834270535828667, loss = 50582.71046448686
21. b = 0.060830413237732354, c = 0.11677514835658626, loss = 23836.651815294274
22. b = 0.05196789447732566, c = 0.11507068215568887, loss = 3723.6422274577717
23. b = 0.05162219593181882, c = 0.11261907551667362, loss = 3632.4460192722363
24. b = 0.04857207270227998, c = 0.11283116783827159, loss = 5821.098745018263
25. b = 0.053134632876166016, c = 0.11242893153633476, loss = 4431.646369650236
26. b = 0.05085804720915753, c = 0.11265381755660105, loss = 3679.943828584574
27. b = 0.05160314914411082, c = 0.11241079505360786, loss = 3629.787082033925
28. b = 0.052347127340599924, c = 0.11226599714163608, loss = 3870.2546740432335
29. b = 0.051975138242355376, c = 0.11233839609762197, loss = 3714.5084536922413
30. b = 0.05122061961067981, c = 0.1124194949597824, loss = 3618.8228191219728

Результаты оценки:
b = 0.05122061961067981
c = 0.1124194949597824

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

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

Result:

In 30 iterations, we have selected the following parameters:

print("Evaluation results: b = ", b, ", c = ", c)

Результаты оценки: b = 0.05122061961067981, c = 0.1124194949597824

We are quite satisfied with the results.

To improve the result, you can experiment with initial assumptions, optimization settings, and submit various input data to the experimental bench.

Row	time	Ступенчатая функция	Add
	Float64	Float64	Float64
1	0.0	0.0	0.197934
2	0.01	0.0	0.631221
3	0.02	0.0	-0.205922
4	0.03	0.0	-1.02133
5	0.04	0.0	0.48395
6	0.05	0.0	-0.200081
7	0.06	0.0	-0.234782
8	0.07	0.0	-0.52721
9	0.08	0.0	-0.632835
10	0.09	0.0	0.791541
11	0.1	0.0	-0.18815
12	0.11	0.0	-0.39816
13	0.12	0.0	-0.066747
⋮	⋮	⋮	⋮
990	9.89	12.0	50.9197
991	9.9	12.0	53.021
992	9.91	12.0	51.6726
993	9.92	12.0	51.0411
994	9.93	12.0	51.6731
995	9.94	12.0	51.402
996	9.95	12.0	51.2705
997	9.96	12.0	52.398
998	9.97	12.0	51.0498
999	9.98	12.0	52.5656
1000	9.99	12.0	52.7554
1001	10.0	12.0	51.5298