updates

Author: baggepinnen

ext/ControlSystemIdentificationLSOptExt.jl

@@ -178,7 +178,13 @@ function _inner_pem(
         resid = zeros(T * ny)
         J = ForwardDiff.jacobian(residuals!, resid, res.minimizer)
         residuals!(resid, res.minimizer)
-        (T - length(p_guess)) / dot(resid, resid) * Symmetric(J' * J)
+        Σout = vec(sum(abs2, reshape(resid, ny, T), dims=2)) ./ (T - length(p_guess))
+        inds = range(1, step=ny, length=T)
+        for i = 1:ny
+            @views J[inds, :] ./= sqrt(Σout[i])
+            inds = inds .+ 1
+        end
+        Symmetric(J' * J)
     end
 
     ukf = get_ukf(res.minimizer)

src/pem.jl

@@ -642,10 +642,10 @@ using ChainRules
 using ForwardDiffChainRules
 @ForwardDiff_frule LinearAlgebra.exp!(x1::AbstractMatrix{<:ForwardDiff.Dual}) true
 
-function inner_constructor(p, constructor, Ts)
+function inner_constructor(p, constructor, Ts, method)
     sys = constructor(p)
     if iscontinuous(sys)
-        return c2d(sys, Ts, :zoh), p.K, p.x0
+        return c2d(sys, Ts, method), p.K, p.x0
     end
     return sys, p.K, p.x0
 end
@@ -710,6 +710,7 @@ function structured_pem(
         # alphaguess = LineSearches.InitialStatic(alpha = 0.95),
         linesearch = LineSearches.BackTracking(),
     ),
+    method = :zoh,
     store_trace = true,
     show_trace = true,
     show_every = 50,
@@ -742,22 +743,34 @@ function structured_pem(
     end
     p0ca = ComponentArray(p = p0, K = T.(K0), x0 = x0i)
     function predloss(p)
-        sysi, Ki, x0 = inner_constructor(p, constructor, d.Ts)
+        sysi, Ki, x0 = inner_constructor(p, constructor, d.Ts, method)
         syso = PredictionStateSpace(sysi, Ki, 0, 0)
         Pyh = ControlSystemsBase.observer_predictor(syso; h)
         simres = lsim(Pyh, pd; x0=Vector(x0))
-        yh = simres.y 
+        yh = simres.y
         yh .= metric.(yh .- d.y)
         c1 = sum(yh)
         c1 + inner_regularizer(regularizer, p, syso, simres)
     end
     function simloss(p)
-        syssim, _, x0 = inner_constructor(p, constructor, d.Ts)
+        syssim, _, x0 = inner_constructor(p, constructor, d.Ts, method)
         simres = lsim(syssim, d; x0)
-        y = simres.y 
+        y = simres.y
         y .= metric.(y .- d.y)
         sum(y) + inner_regularizer(regularizer, p, syssim, simres)
     end
+    function simresiduals(p)
+        syssim, _, x0 = inner_constructor(p, constructor, d.Ts, method)
+        simres = lsim(syssim, d; x0)
+        vec(simres.y .- d.y)
+    end
+    function predresiduals(p)
+        sysi, Ki, x0 = inner_constructor(p, constructor, d.Ts, method)
+        syso = PredictionStateSpace(sysi, Ki, 0, 0)
+        Pyh = ControlSystemsBase.observer_predictor(syso; h)
+        simres = lsim(Pyh, pd; x0=Vector(x0))
+        vec(simres.y .- d.y[:, 1:size(simres.y, 2)])
+    end
     res = Optim.optimize(
         pred ? predloss : simloss,
         p0ca,
@@ -767,7 +780,7 @@ function structured_pem(
             time_limit, x_abstol, f_abstol, g_tol, f_calls_limit, g_calls_limit);
         autodiff = AutoForwardDiff(),
     )
-    sys_opt, K_opt, x0_opt = inner_constructor(res.minimizer, constructor, d.Ts)
+    sys_opt, K_opt, x0_opt = inner_constructor(res.minimizer, constructor, d.Ts, method)
     sysp_opt = PredictionStateSpace(
         sys_opt,
         pred ? K_opt : zeros(T, nx, ny),
@@ -776,5 +789,20 @@ function structured_pem(
         zeros(T, nx, ny),
     )
     pred && !isstable(observer_predictor(sysp_opt)) && @warn("Estimated predictor dynamics A-KC is unstable")
-    (; sys=sysp_opt, x0=x0_opt, res)
+    function Λ()
+        residuals = pred ? predresiduals : simresiduals
+        resid = residuals(res.minimizer)
+        J = ForwardDiff.jacobian(residuals, res.minimizer)
+        Tn = length(resid) ÷ ny
+        npar = length(res.minimizer)
+        Σout = vec(sum(abs2, reshape(resid, ny, Tn), dims=2)) ./ (Tn - npar)
+        inds = range(1, step=ny, length=Tn)
+        for i = 1:ny
+            @views J[inds, :] ./= sqrt(Σout[i])
+            inds = inds .+ 1
+        end
+        Symmetric(J' * J)
+    end
+    ax = ComponentArrays.getaxes(res.minimizer)
+    (; sys=sysp_opt, x0=x0_opt, res, Λ, constructor=p->inner_constructor(ComponentArray(p, ax), constructor, d.Ts, method))
 end

test/test_nonlinear_pem.jl

@@ -2,6 +2,7 @@ using ControlSystemIdentification, SeeToDee, LowLevelParticleFilters
 using LeastSquaresOptim
 using StaticArrays
 using LinearAlgebra
+using Random, Statistics
 using Test
 
 function quadtank(h, u, p, t)
@@ -110,4 +111,73 @@ modelγ = ControlSystemIdentification.nonlinear_pem(d, discrete_dynamics, measur
 
 γ = 1e4*ones(8) # Strong regularization to more or less force the optimum to be at the initial guess
 modelγ = ControlSystemIdentification.nonlinear_pem(d, discrete_dynamics, measurement, p0, x00, R1, 100R2, nu; γ, iterations=300, show_every=50, optimizer=Dogleg())
-@test norm(modelγ.p - p0) < 0.01norm(model.p - p0)
+@test norm(modelγ.p - p0) < 0.01norm(model.p - p0)
+
+
+## Covariance scaling vs bootstrap (ny=1 and ny=3)
+
+@testset "nonlinear_pem covariance scaling" begin
+    @info "Testing nonlinear_pem covariance scaling against bootstrap (ny=1 and ny=3)"
+
+    Ts_b = 1.0
+    nT_b = 200
+    σ_b = 0.3
+    N_reps = 30
+
+    dyn_bs(x, u, p, t) = SA[p[1] * x[1] + p[2] * u[1]]
+    meas_1(x, u, p, t) = SA[x[1]]
+    meas_3(x, u, p, t) = SA[x[1], 0.5 * x[1], -0.3 * x[1]]
+
+    p_true_bs = [0.9, 1.0]
+    x0_bs = SA[1.0]
+
+    Random.seed!(0)
+    u_vec = [SA[0.5 * randn()] for _ in 1:nT_b]
+    U_bs = reduce(hcat, u_vec)
+
+    function rollout_manual(f, x0, uvec, p)
+        x = x0
+        traj = Vector{typeof(x0)}(undef, length(uvec))
+        for k in eachindex(uvec)
+            traj[k] = x
+            x = f(x, uvec[k], p, 0)
+        end
+        traj
+    end
+    x_traj = rollout_manual(dyn_bs, x0_bs, u_vec, p_true_bs)
+
+    function bootstrap_compare(meas, ny, seed_offset)
+        y_clean_vec = meas.(x_traj, u_vec, Ref(p_true_bs), 0)
+        Y_clean = reduce(hcat, y_clean_vec)
+
+        R1_bs = Matrix(1e-6 * I, 1, 1)
+        R2_bs = Matrix(σ_b^2 * I, ny, ny)
+
+        p_ests   = zeros(2, N_reps)
+        var_ests = zeros(2, N_reps)
+
+        redirect_stdout(devnull) do
+            for rep in 1:N_reps
+                Random.seed!(seed_offset + rep)
+                Y_noisy = Y_clean .+ σ_b .* randn(ny, nT_b)
+                d_rep = iddata(Y_noisy, U_bs, Ts_b)
+                model_b = ControlSystemIdentification.nonlinear_pem(
+                    d_rep, dyn_bs, meas, p_true_bs, [x0_bs...], R1_bs, R2_bs, 1,
+                )
+                p_ests[:, rep]   = model_b.p
+                var_ests[:, rep] = diag(inv(model_b.Λ()))[1:2]
+            end
+        end
+        vec(mean(var_ests, dims=2)), vec(var(p_ests, dims=2))
+    end
+
+    # ny = 1
+    mean_est_1, emp_var_1 = bootstrap_compare(meas_1, 1, 100)
+    @test 0.6 < mean_est_1[1] / emp_var_1[1] < 1.7
+    @test 0.6 < mean_est_1[2] / emp_var_1[2] < 1.7
+
+    # ny = 3 — catches the (T·ny vs T) bug in the previous master formula
+    mean_est_3, emp_var_3 = bootstrap_compare(meas_3, 3, 200)
+    @test 0.6 < mean_est_3[1] / emp_var_3[1] < 1.7
+    @test 0.6 < mean_est_3[2] / emp_var_3[2] < 1.7
+end

test/test_pem.jl

@@ -1,6 +1,6 @@
 using ControlSystemIdentification, Optim, ControlSystemsBase
 using ControlSystemsBase.DemoSystems: resonant
-using Random, Test
+using Random, Test, Statistics, LinearAlgebra
 Random.seed!(1)
 T = 1000
 sys = c2d(resonant(ω0 = 0.1) * tf(1, [0.1, 1]), 1)# generate_system(nx, nu, ny)
@@ -358,7 +358,57 @@ end
     res = ControlSystemIdentification.structured_pem(d, nx; focus=:simulation, p0, constructor, regularizer, show_trace = false, show_every = 1, iterations=300)
     p_est = res.res.minimizer.p
     @test norm(p_est - p_true) < 1e-6
-    
+
+end
+
+
+@testset "structured_pem covariance scaling" begin
+    @info "Testing structured_pem covariance scaling against bootstrap (ny=1 and ny=3)"
+
+    Ts_b = 0.1
+    nT_b = 200
+    t_b = 0:Ts_b:Ts_b*(nT_b-1)
+    σ_b = 0.3
+    N_reps = 50
+
+    Random.seed!(0)
+    u_b = randn(1, nT_b)
+
+    function bootstrap_compare(constructor, p_true, seed_offset)
+        y_clean = lsim(c2d(constructor(p_true), Ts_b), u_b, t_b).y
+        p_ests   = zeros(length(p_true), N_reps)
+        var_ests = zeros(length(p_true), N_reps)
+        for rep in 1:N_reps
+            Random.seed!(seed_offset + rep)
+            d_rep = iddata(y_clean .+ σ_b .* randn(size(y_clean)), u_b, Ts_b)
+            res = ControlSystemIdentification.structured_pem(
+                d_rep, 1; focus=:simulation, p0=p_true, constructor,
+                show_trace=false, iterations=100,
+            )
+            p_ests[:, rep]   = collect(res.res.minimizer.p)
+            var_ests[:, rep] = diag(inv(res.Λ()))[1:length(p_true)]
+        end
+        vec(mean(var_ests, dims=2)), vec(var(p_ests, dims=2))
+    end
+
+    # --- ny = 1 (catches the missing factor of ½ — buggy formula would give ratio ≈ 0.5) ---
+    constructor_1 = function (p)
+        a, b = p
+        ss(reshape([a], 1, 1), reshape([b], 1, 1), reshape([1.0], 1, 1), 0)
+    end
+    mean_est_1, emp_var_1 = bootstrap_compare(constructor_1, [-1.0, 1.0], 100)
+    @test 0.6 < mean_est_1[1] / emp_var_1[1] < 1.7
+    @test 0.6 < mean_est_1[2] / emp_var_1[2] < 1.7
+
+    # --- ny = 3 (catches the T vs T·ny error — buggy formula would give ratio ≈ 1.5) ---
+    C_3 = reshape([1.0, 0.5, -0.3], 3, 1)
+    constructor_3 = function (p)
+        a, b = p
+        ss(reshape([a], 1, 1), reshape([b], 1, 1), C_3, zeros(3, 1))
+    end
+    mean_est_3, emp_var_3 = bootstrap_compare(constructor_3, [-1.0, 1.0], 200)
+    @test 0.6 < mean_est_3[1] / emp_var_3[1] < 1.7
+    @test 0.6 < mean_est_3[2] / emp_var_3[2] < 1.7
 end