Школа системного моделирования

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

EngeePkg.purge()

Загрузка данных из таблицы

using DataFrames, CSV;

Import_data = "/user/Школа моделирования/Лекция 12/exp.csv" # @param {type:"filepicker"}
table = CSV.read(Import_data, DataFrame);
print(".csv файл импортирован")
Показать_таблицу = 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, Столбец_выхода];

df_in = DataFrame(time=exp_time_data, value=exp_in_data)
df_in = Float64.(df_in);
wa_in = WorkspaceArray(string(rand()), df_in)

WorkspaceArray{Float64}("0.4774595750319828")

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

input = plot(exp_time_data, exp_in_data, title = "Вход", xlabel="t, c",legend=false)
output_measure = plot(exp_time_data, exp_out_data, title = "Выход", xlabel="t, c", legend=false)
plot(output_measure,input,layout = (2, 1))

Запуск модели, в которой оцениваем параметры

model_name = "parameter_estimation" # @param {type:"string"}
if model_name in [m.name for m in engee.get_all_models()]
    m = engee.open( model_name ) # загрузка модели
else
    m = engee.load( model_name * ".engee" )
end
logout = engee.run(m, verbose=false)
output = plot(exp_time_data, exp_out_data, title = "Выход", xlabel="t, c", label = "эксперимент", legend=true)
plot!(logout["w"].time, logout["w"].value, label = "моделирование", linewidth=2)
plot(output,input,layout = (2, 1))

Подключение библиотек для оценки

import Pkg

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

using DataFrames, NLopt, Interpolations

[ Info: Precompiling NLopt [76087f3c-5699-56af-9a33-bf431cd00edd] (cache misses: wrong julia version (2))
GKS: cannot open display - headless operation mode active
GKS: cannot open display - headless operation mode active

# @markdown # Оценка параметров модели
# @markdown ## Для начала оценки параметров модели введите:
# @markdown ### Название модели:
# Открываем модель
model_name = "parameter_estimation" # @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 = [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_name = "w" # @param {type:"string"}
output_name in keys(engee.run(model, verbose=false)) || throw(ArgumentError("Сигнал $output_name не найден среди логируемых сигналов."))

# @markdown ### Значение выходного сигнала:
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 ### Имена оцениваемых параметров:
# Создаем вектор переменных, их начальные значения и ограничения
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 ### Начальные значения параметров:
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 ### Нижняя граница значений параметров:
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 ### Верхняя граница значений параметров:
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 ## Настройки оптимизации:
# @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. b = 0.01, c = 0.7, loss = 780498.8498020292
2. b = 0.7575756087138917, c = 3.2108789930235626, loss = 2.1245193607542957e6
3. b = 0.25857951120956724, c = 5.352505475679727, loss = 2.0738436200875745e6
4. b = 0.4879434824253316, c = 2.401502081602793, loss = 1.9092307050506342e6
5. b = 0.5701990541362018, c = 1.4620962731302711, loss = 1.8338366245770543e6
6. b = 0.01, c = 0.7, loss = 780498.8498020292
7. b = 0.0175, c = 0.7, loss = 332787.244494691
8. b = 0.0175, c = 1.2249999999999999, loss = 95303.07549796053
9. b = 0.024125589953757652, c = 1.4710134814331561, loss = 372.54975185700516
10. b = 0.0235427244628947, c = 1.9944256634807744, loss = 47688.155289581984
11. b = 0.02772267541502742, c = 1.3968175697993195, loss = 1835.6069863879695
12. b = 0.023494942153417096, c = 1.3474102819420608, loss = 7106.005910872096
13. b = 0.02445972994413676, c = 1.6001625417760958, loss = 3401.975962986597
14. b = 0.024292659948947207, c = 1.5355880116046259, loss = 777.3210473849299
15. b = 0.023873666305769287, c = 1.47332582857537, loss = 485.93642217047756
16. b = 0.02469023141460299, c = 1.4186262432314458, loss = 756.5449100133457
17. b = 0.023670853041123006, c = 1.4630504152109778, loss = 863.908298359532
18. b = 0.024262137263233792, c = 1.4843477874596862, loss = 265.83891636875177
19. b = 0.024490650781413815, c = 1.4807015919209716, loss = 281.71426092223356
20. b = 0.024284405976300036, c = 1.5006798158307284, loss = 306.46997086063635
21. b = 0.02419815504304029, c = 1.4774752326427887, loss = 297.7354537901099
22. b = 0.024262862611504434, c = 1.4884490356725903, loss = 263.6365894081367
23. b = 0.024321308127884, c = 1.4887406018838092, loss = 270.4368125814954
24. b = 0.024207185565888086, c = 1.4897269798182549, loss = 260.27388471112744
25. b = 0.024161375816942117, c = 1.4922843227281501, loss = 257.4601841499203
26. b = 0.02410282120019854, c = 1.4924342011644545, loss = 256.2355731036028
27. b = 0.024083154949782214, c = 1.496297837282771, loss = 252.78322413705152
28. b = 0.024045623526730214, c = 1.5040688785228957, loss = 256.73213071926546
29. b = 0.02403402025138868, c = 1.4940633029086987, loss = 255.46987339148768
30. b = 0.02408347083677361, c = 1.4983484993199812, loss = 253.31001020607067

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

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

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

Row	time	Ступенчатая функция	Физическое моделирование
	Float64	Float64	Float64
1	0.0	0.0	0.0
2	0.01	0.0	0.0
3	0.02	0.0	0.0
4	0.03	0.0	0.0
5	0.04	0.0	0.0
6	0.05	0.0	0.0
7	0.06	0.0	0.0
8	0.07	0.0	0.0
9	0.08	0.0	0.0
10	0.09	0.0	0.0
11	0.1	0.0	0.0
12	0.11	0.0	0.0
13	0.12	0.0	0.0
⋮	⋮	⋮	⋮
990	9.89	12.0	51.6667
991	9.9	12.0	51.6667
992	9.91	12.0	51.6667
993	9.92	12.0	51.6667
994	9.93	12.0	51.6667
995	9.94	12.0	51.6667
996	9.95	12.0	51.6667
997	9.96	12.0	51.6667
998	9.97	12.0	51.6667
999	9.98	12.0	51.6667
1000	9.99	12.0	51.6667
1001	10.0	12.0	51.6667

Модуль 4. Лекция 12 - Школа системного моделирования

Школа системного моделирования

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

Загрузка данных из таблицы

Запуск модели, в которой оцениваем параметры

Подключение библиотек для оценки

Type

Short description

Related posts

Tags

Categories

Languages

Formats

Level

Status

Published

Updated

Source