Model adaptation and tuning online changes

[caxelrud]

Hi. I tested the model adaptation and tuning online changes using the code that follows my comments. I am using:

M_Hp1=diagm(repeat(Mwt,Hp))
Ñ_Hc1=diagm(push!(repeat(Nwt,Hc),1e5))

Let me know if it correct.

One suggestion, to make it simpler, is to use Mwt and Nwt parameters as in LinMPC:

mpc = LinMPC(estim, Hp=200, Hc=50, Mwt=[1.0, 1.0, 0.01], Nwt=[2.0, 2.0]) Regards, Celso

using ModelPredictiveControl, ControlSystemsBase
using Plots
using LinearAlgebra

Hp=200;Hc=50
Ny=3;Nu=2;Nd=1
G = [ tf(1.90, [180, 1]) tf(1.90, [180,1]) tf(1.90, [180,1]);
      tf(-0.74,[80, 1])  tf(0.74, [80, 1]) tf(0.74, [80, 1]);
      tf(-0.5,[100, 1])  tf(2.0, [250, 1]) tf(0.1, [80, 1])]
Ts = 10.0
model = setop!(LinModel(G, Ts,i_u=[1,2],i_d=[3]), uop=[20, 20], yop=[50, 30, 25], dop=[20])
setname!(model; u=["mv1","mv2"], y=["cv1","cv2","ccv1"],d=["dv1"])
estim = KalmanFilter(model)
mpc = LinMPC(estim, Hp=200, Hc=50, Mwt=[1.0, 1.0, 0.01], Nwt=[2.0, 2.0])
umin=[10,10];umax=[55,55];Δumax=[0.5,0.5]
ymin=[20,20,20];ymax=[70,70,35]
mpc = setconstraint!(mpc, umin=umin,umax=umax,Δumax=Δumax,c_ymin=fill(0.0,Ny),c_ymax=fill(0.0,Ny))
mpc = setconstraint!(mpc, ymin=ymin,ymax=ymax) 

G1 = [ tf(3.00, [180, 1]) tf(1.90, [180,1]) tf(1.90, [180,1]);
      tf(-0.74,[80, 1])  tf(0.74, [80, 1]) tf(0.74, [80, 1]);
      tf(-0.5,[100, 1])  tf(2.0, [250, 1]) tf(0.1, [80, 1])]
model1 = setop!(LinModel(G1, Ts,i_u=[1,2],i_d=[3]), uop=[20, 20], yop=[50, 30, 25], dop=[20])

Mwt=[1.0, 1.0, 0.01]; Nwt=[1.0, 1.0]
M_Hp1=diagm(repeat(Mwt,Hp))
Ñ_Hc1=diagm(push!(repeat(Nwt,Hc),1e5))

u=model.uop; d=model.dop; y = model(d)  # x->y 
x=initstate!(mpc,u,y,d) # u,y,d->x
ry, d = [50.0,30.0,25.0], [20]
y = model(d) 
u = mpc(ry,d) 

function test_mpc(mpc, model)
    N = 200
    ry, d = [50.0,30.0,25.0], [20]
    u_data, y_data, ry_data = zeros(model.nu, N), zeros(model.ny, N), zeros(model.ny, N)
    for i = 1:N
        i == 51  && (ry = [50, 40, 25])
        i == 101 && (ry = [54, 30, 25] ;setmodel!(mpc, model1;M_Hp=M_Hp1,Ñ_Hc=Ñ_Hc1))
        i == 151 && (d = [25]) 
        y = model(d) 
        u = mpc(ry,d) 
        u_data[:,i], y_data[:,i], ry_data[:,i] = u, y, ry
        updatestate!(mpc, u, y, d) # update mpc state estimate
        updatestate!(model, u , d) # update simulator with the load disturbance
    end
    return u_data, y_data, ry_data
end

function plot_data(t_data, u_data, y_data, ry_data)
    p1 = plot(t_data, y_data[1,:], label="meas.", ylabel="level")
    plot!(p1, t_data, ry_data[1,:], label="setpoint", linestyle=:dash, linetype=:steppost)
    plot!(p1, t_data, fill(45,size(t_data)), label="min", linestyle=:dot, linewidth=1.5)
    hline!(p1,[ymin[1],ymax[1]],linestyle=:dot,label="")

    p2 = plot(t_data, y_data[2,:], label="meas.", legend=:topleft, ylabel="temp.")
    plot!(p2, t_data, ry_data[2,:],label="setpoint", linestyle=:dash, linetype=:steppost)
    hline!(p2,[ymin[2],ymax[2]],linestyle=:dot,label="")

    p3 = plot(t_data, y_data[3,:], label="meas.", legend=:topleft, ylabel="temp.2")
    plot!(p3, t_data, ry_data[3,:],label="setpoint", linestyle=:dash, linetype=:steppost)
    hline!(p3,[ymin[3],ymax[3]],linestyle=:dot,label="")

    p4 = plot(t_data,u_data[1,:],label="cold", linetype=:steppost, ylabel="flow rate")
    plot!(p4, t_data,u_data[2,:],label="hot", linetype=:steppost, xlabel="time (s)")
    hline!(p4,[umin[1],umax[1]],linestyle=:dot,label="")

    return plot(p1, p2, p3, p4, layout=(4,1),size=(900,500))
end

u_data, y_data, ry_data = test_mpc(mpc, model)
t_data = Ts*(0:(size(y_data,2)-1))
plot_data(t_data, u_data, y_data, ry_data)

[franckgaga]

yes good point that would simplify the interface!

I'll take a look at your code soon.

[franckgaga]

Your code seems ok to me. Some minor comments : you should avoid global variable if possible (for performance and readability). Here model1, M_Hp1 and Ñ_Hc1 are global variables in test_mpc function. Add new arguments in test_mpc function to avoid globals.

Also, I don't know if it was intended but, at i==101, you modified the plant model in the controller but not the plant simulator. About that, be careful if you assign a new plant simulator during simulations. The simulators states are stored inside the objects (at obj.x0). Thus, at i==101, if you do model=model1, the state of the plant simulator will jump to model1 internal value at this moment, that is, zero.

I'll work on your suggestion to simplify the interface.