Skip to content

Commit ebf2a42

Browse files
franckgagafranckgaga
authored
Merge pull request #222 from JuliaControl/new_bench
bench: adding some `PredictiveController` benchmarks
2 parents 69cf0bf + bd3a48e commit ebf2a42

File tree

7 files changed

+827
-138
lines changed

7 files changed

+827
-138
lines changed

‎.gitignore‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,3 +9,4 @@ lcov.info
99
benchmark/Manifest.toml
1010
/.benchmarkci
1111
/benchmark/*.json
12+
results_ModelPredictiveControl*

‎benchmark/0_bench_setup.jl‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
Ts = 400.0
2+
sys = [ tf(1.90,[1800.0,1]) tf(1.90,[1800.0,1]) tf(1.90,[1800.0,1]);
3+
tf(-0.74,[800.0,1]) tf(0.74,[800.0,1]) tf(-0.74,[800.0,1]) ]
4+
function f_lin!(ẋ, x, u, d, p)
5+
mul!(ẋ, p.A, x)
6+
mul!(ẋ, p.Bu, u, 1, 1)
7+
mul!(ẋ, p.Bd, d, 1, 1)
8+
return nothing
9+
end
10+
function h_lin!(y, x, d, p)
11+
mul!(y, p.C, x)
12+
mul!(y, p.Dd, d, 1, 1)
13+
return nothing
14+
end
15+
linmodel = setop!(LinModel(sys, Ts, i_d=[3]), uop=[10, 50], yop=[50, 30], dop=[5])
16+
nonlinmodel = NonLinModel(f_lin!, h_lin!, Ts, 2, 4, 2, 1, p=linmodel, solver=nothing)
17+
nonlinmodel = setop!(nonlinmodel, uop=[10, 50], yop=[50, 30], dop=[5])
18+
u, d, y = [10, 50], [5], [50, 30]

‎benchmark/1_bench_sim_model.jl‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
## ----------------- Unit tests ----------------------------------------------------------
2+
const UNIT_MODEL = SUITE["unit tests"]["SimModel"]
3+
4+
UNIT_MODEL["LinModel"]["updatestate!"] =
5+
@benchmarkable(
6+
updatestate!($linmodel, $u, $d);
7+
)
8+
UNIT_MODEL["LinModel"]["evaloutput"] =
9+
@benchmarkable(
10+
evaloutput($linmodel, $d);
11+
)
12+
UNIT_MODEL["NonLinModel"]["updatestate!"] =
13+
@benchmarkable(
14+
updatestate!($nonlinmodel, $u, $d);
15+
)
16+
UNIT_MODEL["NonLinModel"]["evaloutput"] =
17+
@benchmarkable(
18+
evaloutput($nonlinmodel, $d);
19+
)
20+
UNIT_MODEL["NonLinModel"]["linearize!"] =
21+
@benchmarkable(
22+
linearize!($linmodel, $nonlinmodel);
23+
)
24+
25+
## ----------------- Case studies ---------------------------------------------------------
26+
const CASE_MODEL = SUITE["case studies"]["SimModel"]
27+
# TODO: Add case study benchmarks for SimModel

‎benchmark/2_bench_state_estim.jl‎

Lines changed: 285 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,285 @@
1+
## ----------------- Unit tests -----------------------------------------------------------
2+
const UNIT_ESTIM = SUITE["unit tests"]["StateEstimator"]
3+
4+
skf = SteadyKalmanFilter(linmodel)
5+
UNIT_ESTIM["SteadyKalmanFilter"]["preparestate!"] =
6+
@benchmarkable(
7+
preparestate!($skf, $y, $d),
8+
)
9+
UNIT_ESTIM["SteadyKalmanFilter"]["updatestate!"] =
10+
@benchmarkable(
11+
updatestate!($skf, $u, $y, $d),
12+
setup=preparestate!($skf, $y, $d),
13+
)
14+
UNIT_ESTIM["SteadyKalmanFilter"]["evaloutput"] =
15+
@benchmarkable(
16+
evaloutput($skf, $d),
17+
setup=preparestate!($skf, $y, $d),
18+
)
19+
20+
kf = KalmanFilter(linmodel, nint_u=[1, 1], direct=false)
21+
UNIT_ESTIM["KalmanFilter"]["preparestate!"] =
22+
@benchmarkable(
23+
preparestate!($kf, $y, $d),
24+
)
25+
UNIT_ESTIM["KalmanFilter"]["updatestate!"] =
26+
@benchmarkable(
27+
updatestate!($kf, $u, $y, $d),
28+
setup=preparestate!($kf, $y, $d),
29+
)
30+
UNIT_ESTIM["KalmanFilter"]["evaloutput"] =
31+
@benchmarkable(
32+
evaloutput($kf, $d),
33+
setup=preparestate!($kf, $y, $d),
34+
)
35+
36+
lo = Luenberger(linmodel, nint_u=[1, 1])
37+
UNIT_ESTIM["Luenberger"]["preparestate!"] =
38+
@benchmarkable(
39+
preparestate!($lo, $y, $d),
40+
)
41+
UNIT_ESTIM["Luenberger"]["updatestate!"] =
42+
@benchmarkable(
43+
updatestate!($lo, $u, $y, $d),
44+
setup=preparestate!($lo, $y, $d),
45+
)
46+
UNIT_ESTIM["Luenberger"]["evaloutput"] =
47+
@benchmarkable(
48+
evaloutput($lo, $d),
49+
setup=preparestate!($lo, $y, $d),
50+
)
51+
52+
im_lin = InternalModel(linmodel)
53+
im_nonlin = InternalModel(nonlinmodel)
54+
UNIT_ESTIM["InternalModel"]["preparestate!"]["LinModel"] =
55+
@benchmarkable(
56+
preparestate!($im_lin, $y, $d),
57+
)
58+
UNIT_ESTIM["InternalModel"]["updatestate!"]["LinModel"] =
59+
@benchmarkable(
60+
updatestate!($im_lin, $u, $y, $d),
61+
setup=preparestate!($im_lin, $y, $d),
62+
)
63+
UNIT_ESTIM["InternalModel"]["evaloutput"]["LinModel"] =
64+
@benchmarkable(
65+
evaloutput($im_lin, $d),
66+
setup=preparestate!($im_lin, $y, $d),
67+
)
68+
UNIT_ESTIM["InternalModel"]["preparestate!"]["NonLinModel"] =
69+
@benchmarkable(
70+
preparestate!($im_nonlin, $y, $d),
71+
)
72+
UNIT_ESTIM["InternalModel"]["updatestate!"]["NonLinModel"] =
73+
@benchmarkable(
74+
updatestate!($im_nonlin, $u, $y, $d),
75+
setup=preparestate!($im_nonlin, $y, $d),
76+
)
77+
UNIT_ESTIM["InternalModel"]["evaloutput"]["NonLinModel"] =
78+
@benchmarkable(
79+
evaloutput($im_nonlin, $d),
80+
setup=preparestate!($im_nonlin, $y, $d),
81+
)
82+
83+
ukf_lin = UnscentedKalmanFilter(linmodel)
84+
ukf_nonlin = UnscentedKalmanFilter(nonlinmodel)
85+
UNIT_ESTIM["UnscentedKalmanFilter"]["preparestate!"]["LinModel"] =
86+
@benchmarkable(
87+
preparestate!($ukf_lin, $y, $d),
88+
)
89+
UNIT_ESTIM["UnscentedKalmanFilter"]["updatestate!"]["LinModel"] =
90+
@benchmarkable(
91+
updatestate!($ukf_lin, $u, $y, $d),
92+
setup=preparestate!($ukf_lin, $y, $d),
93+
)
94+
UNIT_ESTIM["UnscentedKalmanFilter"]["evaloutput"]["LinModel"] =
95+
@benchmarkable(
96+
evaloutput($ukf_lin, $d),
97+
setup=preparestate!($ukf_lin, $y, $d),
98+
)
99+
UNIT_ESTIM["UnscentedKalmanFilter"]["preparestate!"]["NonLinModel"] =
100+
@benchmarkable(
101+
preparestate!($ukf_nonlin, $y, $d),
102+
)
103+
UNIT_ESTIM["UnscentedKalmanFilter"]["updatestate!"]["NonLinModel"] =
104+
@benchmarkable(
105+
updatestate!($ukf_nonlin, $u, $y, $d),
106+
setup=preparestate!($ukf_nonlin, $y, $d),
107+
)
108+
UNIT_ESTIM["UnscentedKalmanFilter"]["evaloutput"]["NonLinModel"] =
109+
@benchmarkable(
110+
evaloutput($ukf_nonlin, $d),
111+
setup=preparestate!($ukf_nonlin, $y, $d),
112+
)
113+
114+
ekf_lin = ExtendedKalmanFilter(linmodel, nint_u=[1, 1], direct=false)
115+
ekf_nonlin = ExtendedKalmanFilter(nonlinmodel, nint_u=[1, 1], direct=false)
116+
UNIT_ESTIM["ExtendedKalmanFilter"]["preparestate!"]["LinModel"] =
117+
@benchmarkable(
118+
preparestate!($ekf_lin, $y, $d),
119+
)
120+
UNIT_ESTIM["ExtendedKalmanFilter"]["updatestate!"]["LinModel"] =
121+
@benchmarkable(
122+
updatestate!($ekf_lin, $u, $y, $d),
123+
setup=preparestate!($ekf_lin, $y, $d),
124+
)
125+
UNIT_ESTIM["ExtendedKalmanFilter"]["evaloutput"]["LinModel"] =
126+
@benchmarkable(
127+
evaloutput($ekf_lin, $d),
128+
setup=preparestate!($ekf_lin, $y, $d),
129+
)
130+
UNIT_ESTIM["ExtendedKalmanFilter"]["preparestate!"]["NonLinModel"] =
131+
@benchmarkable(
132+
preparestate!($ekf_nonlin, $y, $d),
133+
)
134+
UNIT_ESTIM["ExtendedKalmanFilter"]["updatestate!"]["NonLinModel"] =
135+
@benchmarkable(
136+
updatestate!($ekf_nonlin, $u, $y, $d),
137+
setup=preparestate!($ekf_nonlin, $y, $d),
138+
)
139+
UNIT_ESTIM["ExtendedKalmanFilter"]["evaloutput"]["NonLinModel"] =
140+
@benchmarkable(
141+
evaloutput($ekf_nonlin, $d),
142+
setup=preparestate!($ekf_nonlin, $y, $d),
143+
)
144+
145+
mhe_lin_curr = MovingHorizonEstimator(linmodel, He=10, direct=true)
146+
mhe_lin_pred = MovingHorizonEstimator(linmodel, He=10, direct=false)
147+
mhe_nonlin_curr = MovingHorizonEstimator(nonlinmodel, He=10, direct=true)
148+
mhe_nonlin_pred = MovingHorizonEstimator(nonlinmodel, He=10, direct=false)
149+
150+
samples, evals, seconds = 5000, 1, 60
151+
UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["LinModel"]["Current form"] =
152+
@benchmarkable(
153+
preparestate!($mhe_lin_curr, $y, $d),
154+
samples=samples, evals=evals, seconds=seconds,
155+
)
156+
UNIT_ESTIM["MovingHorizonEstimator"]["updatestate!"]["LinModel"]["Current form"] =
157+
@benchmarkable(
158+
updatestate!($mhe_lin_curr, $u, $y, $d),
159+
setup=preparestate!($mhe_lin_curr, $y, $d),
160+
samples=samples, evals=evals, seconds=seconds,
161+
)
162+
UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["LinModel"]["Prediction form"] =
163+
@benchmarkable(
164+
preparestate!($mhe_lin_pred, $y, $d),
165+
samples=samples, evals=evals, seconds=seconds,
166+
)
167+
UNIT_ESTIM["MovingHorizonEstimator"]["updatestate!"]["LinModel"]["Prediction form"] =
168+
@benchmarkable(
169+
updatestate!($mhe_lin_pred, $u, $y, $d),
170+
setup=preparestate!($mhe_lin_pred, $y, $d),
171+
samples=samples, evals=evals, seconds=seconds,
172+
)
173+
UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["NonLinModel"]["Current form"] =
174+
@benchmarkable(
175+
preparestate!($mhe_nonlin_curr, $y, $d),
176+
samples=samples, evals=evals, seconds=seconds,
177+
)
178+
UNIT_ESTIM["MovingHorizonEstimator"]["updatestate!"]["NonLinModel"]["Current form"] =
179+
@benchmarkable(
180+
updatestate!($mhe_nonlin_curr, $u, $y, $d),
181+
setup=preparestate!($mhe_nonlin_curr, $y, $d),
182+
samples=samples, evals=evals, seconds=seconds,
183+
)
184+
UNIT_ESTIM["MovingHorizonEstimator"]["preparestate!"]["NonLinModel"]["Prediction form"] =
185+
@benchmarkable(
186+
preparestate!($mhe_nonlin_pred, $y, $d),
187+
samples=samples, evals=evals, seconds=seconds,
188+
)
189+
UNIT_ESTIM["MovingHorizonEstimator"]["updatestate!"]["NonLinModel"]["Prediction form"] =
190+
@benchmarkable(
191+
updatestate!($mhe_nonlin_pred, $u, $y, $d),
192+
setup=preparestate!($mhe_nonlin_pred, $y, $d),
193+
samples=samples, evals=evals, seconds=seconds,
194+
)
195+
196+
## ----------------- Case studies ---------------------------------------------------
197+
const CASE_ESTIM = SUITE["case studies"]["StateEstimator"]
198+
199+
## ----------------- Case study: CSTR -----------------------------------------------------
200+
G = [ tf(1.90, [18, 1]) tf(1.90, [18, 1]);
201+
tf(-0.74,[8, 1]) tf(0.74, [8, 1]) ]
202+
uop, yop = [20, 20], [50, 30]
203+
model = setop!(LinModel(G, 2.0); uop, yop)
204+
plant = setop!(LinModel(G, 2.0); uop, yop)
205+
plant.A[diagind(plant.A)] .-= 0.1 # plant-model mismatch
206+
function test_mhe(mhe, plant)
207+
plant.x0 .= 0; y = plant()
208+
initstate!(mhe, plant.uop, y)
209+
N = 75; u = [20, 20]; ul = 0
210+
U, Y, Ŷ = zeros(2, N), zeros(2, N), zeros(2, N)
211+
for i = 1:N
212+
i == 26 && (u = [15, 25])
213+
i == 51 && (ul = -10)
214+
y = plant()
215+
preparestate!(mhe, y)
216+
ŷ = evaloutput(mhe)
217+
U[:,i], Y[:,i], Ŷ[:,i] = u, y, ŷ
218+
updatestate!(mhe, u, y)
219+
updatestate!(plant, u+[0,ul])
220+
end
221+
return U, Y, Ŷ
222+
end
223+
He = 10; nint_u = [1, 1]; σQint_u = [1, 2]
224+
225+
optim = JuMP.Model(OSQP.Optimizer, add_bridges=false)
226+
direct = true
227+
mhe_cstr_osqp_curr = MovingHorizonEstimator(model; He, nint_u, σQint_u, optim, direct)
228+
mhe_cstr_osqp_curr = setconstraint!(mhe_cstr_osqp_curr, v̂min=[-1, -0.5], v̂max=[+1, +0.5])
229+
JuMP.unset_time_limit_sec(mhe_cstr_osqp_curr.optim)
230+
231+
optim = JuMP.Model(OSQP.Optimizer, add_bridges=false)
232+
direct = false
233+
mhe_cstr_osqp_pred = MovingHorizonEstimator(model; He, nint_u, σQint_u, optim, direct)
234+
mhe_cstr_osqp_pred = setconstraint!(mhe_cstr_osqp_pred, v̂min=[-1, -0.5], v̂max=[+1, +0.5])
235+
JuMP.unset_time_limit_sec(mhe_cstr_osqp_pred.optim)
236+
237+
optim = JuMP.Model(DAQP.Optimizer, add_bridges=false)
238+
direct = true
239+
mhe_cstr_daqp_curr = MovingHorizonEstimator(model; He, nint_u, σQint_u, optim, direct)
240+
mhe_cstr_daqp_curr = setconstraint!(mhe_cstr_daqp_curr, v̂min=[-1, -0.5], v̂max=[+1, +0.5])
241+
JuMP.set_attribute(mhe_cstr_daqp_curr.optim, "eps_prox", 1e-6) # needed to support hessians H≥0
242+
243+
optim = JuMP.Model(DAQP.Optimizer, add_bridges=false)
244+
direct = false
245+
mhe_cstr_daqp_pred = MovingHorizonEstimator(model; He, nint_u, σQint_u, optim, direct)
246+
mhe_cstr_daqp_pred = setconstraint!(mhe_cstr_daqp_pred, v̂min=[-1, -0.5], v̂max=[+1, +0.5])
247+
JuMP.set_attribute(mhe_cstr_daqp_pred.optim, "eps_prox", 1e-6) # needed to support hessians H≥0
248+
249+
optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
250+
direct = true
251+
mhe_cstr_ipopt_curr = MovingHorizonEstimator(model; He, nint_u, σQint_u, optim, direct)
252+
mhe_cstr_ipopt_curr = setconstraint!(mhe_cstr_ipopt_curr, v̂min=[-1, -0.5], v̂max=[+1, +0.5])
253+
JuMP.unset_time_limit_sec(mhe_cstr_ipopt_curr.optim)
254+
255+
optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer,"sb"=>"yes"), add_bridges=false)
256+
direct = false
257+
mhe_cstr_ipopt_pred = MovingHorizonEstimator(model; He, nint_u, σQint_u, optim, direct)
258+
mhe_cstr_ipopt_pred = setconstraint!(mhe_cstr_ipopt_pred, v̂min=[-1, -0.5], v̂max=[+1, +0.5])
259+
JuMP.unset_time_limit_sec(mhe_cstr_ipopt_pred.optim)
260+
261+
samples, evals = 500, 1
262+
CASE_ESTIM["CSTR"]["MovingHorizonEstimator"]["OSQP"]["Current form"] =
263+
@benchmarkable(test_mhe($mhe_cstr_osqp_curr, $plant);
264+
samples=samples, evals=evals
265+
)
266+
CASE_ESTIM["CSTR"]["MovingHorizonEstimator"]["OSQP"]["Prediction form"] =
267+
@benchmarkable(test_mhe($mhe_cstr_osqp_pred, $plant);
268+
samples=samples, evals=evals
269+
)
270+
CASE_ESTIM["CSTR"]["MovingHorizonEstimator"]["DAQP"]["Current form"] =
271+
@benchmarkable(test_mhe($mhe_cstr_daqp_curr, $plant);
272+
samples=samples, evals=evals
273+
)
274+
CASE_ESTIM["CSTR"]["MovingHorizonEstimator"]["DAQP"]["Prediction form"] =
275+
@benchmarkable(test_mhe($mhe_cstr_daqp_pred, $plant);
276+
samples=samples, evals=evals
277+
)
278+
CASE_ESTIM["CSTR"]["MovingHorizonEstimator"]["Ipopt"]["Current form"] =
279+
@benchmarkable(test_mhe($mhe_cstr_ipopt_curr, $plant);
280+
samples=samples, evals=evals
281+
)
282+
CASE_ESTIM["CSTR"]["MovingHorizonEstimator"]["Ipopt"]["Prediction form"] =
283+
@benchmarkable(test_mhe($mhe_cstr_ipopt_pred, $plant);
284+
samples=samples, evals=evals
285+
)

0 commit comments

Comments
 (0)