Help using nonlinear_pem with transfer functions

[Hofsmo]

Hi, I have been testing nonlinear_pem a bit. I know it is marked as experimental, but it seems like it may be quite useful. I managed to make it work calculating the state space model by hand. However, I have already made functions for several transfer functions describing my system. I therefore attempted to use create the state space system from the transfer functions using the function ss provided by the ControlSystem package.

I am able to simulate the system correctly using LowLevelParticleFilters.rollout, however, the function nonlinear_pem fails. I was wondering if there is something fundamentally flawed with my idea and if I should just put in the work and calculate the state space models? If not, I have provided an example of what I attempted, using a first order transfer function:

using SeeToDee, LowLevelParticleFilters
using ControlSystems
using ControlSystemIdentification
using StaticArrays
using LeastSquaresOptim
using Plots
using LinearAlgebra

servo(Tg) = tf([1], [Tg, 1])
Tg = 0.4

function servo_dyn(h, u, p, t)
    Tg = p[1]
    sys0 = ss(servo(Tg));
    sm = @MArray zeros(2)

    for i = 1:size(sys0.A, 1)
        for j = 1:size(sys0.A, 2)
            sm[i] += sys0.A[i,j]*h[j]
        end
        sm[i] += sys0.B[i] *u[1]
    end
    return sm
end

Ts = 0.1
t=range(0; length=500, step=Ts)
discrete_dynamics = SeeToDee.ForwardEuler(servo_dyn, Ts)

u = repeat(randn(50), inner=10)
x = LowLevelParticleFilters.rollout(discrete_dynamics, [0.0, 0.0], u, [Tg])[1:end-1]

function measurement(h, u, p, t)
    Tg = p[1]
    sys0 = ss(servo(Tg));
    SA[sys0.C[1]*h[1]]
end

yw= measurement.(x, 0, Ref([Tg]), 0.0);
yw = [y[1] for y in yw]

y,t = lsim(servo(Tg), u', Array(t))

plot(t, y', label="Transfer function")
plot!(t, yw, label="State space")

R1 = Diagonal([0.1, 0.1])
R2 = 100*Diagonal(0.5^2 * ones(1))
p0 = [Tg]
x0 = [0.0, 0.0]

d = iddata(yw, u, Ts)
model = ControlSystemIdentification.nonlinear_pem(d, discrete_dynamics, measurement, p0, x0, R1, R2, 1 )

[baggepinnen]

Hello 👋

The problem is this line

sm = @MArray zeros(2)

it creates an array of Float64 which is not compatible with the dual numbers used by the AD engine. If you replace it with

sm = @MArray zeros(eltype(h), 2)

it works

[Hofsmo]

Many thanks for the quick response, it works perfectly now :) . I am looking forward to try with more complicated systems.

[baggepinnen]

It looks like you are estimating a linear system still? If so, there is the function structured_pem which is likely more efficient than nonlinear_pem with ForwardEuler discretization.

[Hofsmo]

Thanks for the tip I will test that function. The example was just a part of the system I am investigating. I work with hydro power plants, so there will be a regulator, a servo, and water and turbine dynamics. The water dynamics are nonlinear, and the servo may have backlash. For what I am investigating now a first order transfer function representation of the water and turbine dynamics should be sufficient, and the servo does not seem to have significant backlash so I will change to work with structured_pem