I need help regarding formulating multi-state stochastic MPC in JuMP. In some way my optimization routine is wrong. Any help is much appreciated!
swing MPC
The swing equation is given as,
M\frac{d^2\delta}{dt^2}=P_\mathrm{a}=P_\mathrm{m}-P_\mathrm{e}.
If w=2\pi f is the angular frequency then,
\frac{d\delta}{dt}=2\pi f
\frac{d\omega}{dt}=\frac{P_\mathrm{m}-P_\mathrm{e}}{M}.
The differential equation is given as,
\left(\begin{array}{c}
\dot{\delta}\\
\dot{\omega}
\end{array}\right)=\left[\begin{array}{cc}
0 & 1\\
0 & 0
\end{array}\right]\left(\begin{array}{c}
\delta\\
\omega
\end{array}\right)+\left[\begin{array}{c}
0\\
1
\end{array}\right]P_{\mathrm{m}}+\left[\begin{array}{ccc}
0 & 0 & 0\\
-\frac{1}{M} & \frac{1}{M} & \frac{1}{M}
\end{array}\right]\left(\begin{array}{c}
\begin{array}{c}
P_{\ell}\\
P_{\mathrm{pv}}\\
P_{\mathrm{w}}
\end{array}\end{array}\right).
The difference equation is given by,
\left(\begin{array}{c}
\delta_{k+1}\\
\omega_{k+1}
\end{array}\right)=\left[\begin{array}{cc}
0 & 1\\
0 & 0
\end{array}\right]\left(\begin{array}{c}
\delta_{k}\\
\omega_{k}
\end{array}\right)+\left[\begin{array}{c}
0\\
\frac{1}{M}
\end{array}\right]P_{\mathrm{m,}k}+\left[\begin{array}{ccc}
0 & 0 & 0\\
-\frac{1}{M} & \frac{1}{M} & \frac{1}{M}
\end{array}\right]\left(\begin{array}{c}
\begin{array}{c}
P_{\ell,k}\\
P_{\mathrm{pv,}k}\\
P_{\mathrm{w},k}
\end{array}\end{array}\right)
Where,
x_{k+1}=Ax_{k}+Bu_{k}+Gw_{k}
A =\left[\begin{array}{cc} 0 & 1\\ 0 & 0 \end{array}\right] \\ B =\left[\begin{array}{c} 0\\ \frac{1}{M} \end{array}\right] \\ G =\left[\begin{array}{ccc} 0 & 0 & 0\\ -\frac{1}{M} & \frac{1}{M} & \frac{1}{M} \end{array}\right] \\ x =\left(\begin{array}{c} \delta\\ \omega \end{array}\right) \\ u =P_{\mathrm{m}} \\ w =\left(\begin{array}{c} P_{\ell}\\ P_{\mathrm{pv}}\\ P_{\mathrm{w}} \end{array}\right)
The measuremnent equation is given as,
y_{k}=Cx_{k}
Where,
C=\left[\begin{array}{cc}
0 & 1\end{array}\right]
Cost function
The cost function is given as,
\begin{equation*}
\begin{aligned}
& \underset{u}{\text{min}} & J = \frac{1}{2} \sum_{k=1}^{N_\mathrm{p}} e_\mathrm{k}^T e_\mathrm{k} + u_\mathrm{k-1}^T p u_\mathrm{k-1} \\
& \text{with} \\
& & e_k = y_k - r_k, \\
& & x_{k+1} = Ax_k + Bu_k+Gw_k, \\
& & y_k = Cx_k, \\
& & 0.98*314 \leq x_\mathrm{k}[2] \leq 1.02*314, \\
& & -Inf \leq u_\mathrm{k} \leq Inf, \
\end{aligned}
\end{equation*}
Since we want frequency of grid to be between {(49-51)}{Hz} we should constriant x_k[2]=\omega_k, the angular frequency, between {(0.98*314-1.02*314)}{rad/s}.
The optimized value of u_k=P_{\mathrm{m},k} can take any value. The negative value tells that excess power can be used for power storage in flywheels, inertial bodies, battery, or to a run pump storage plant and positive means a hydro generator should be added to the grid.
r_k is the reference signal for angular frequency which is always 314\ rad/sec.
The parameters and other settings look like this,
using Plots;pyplot();
using LaTeXStrings;
using Optim;
using JuMP;
using Ipopt;
# swing system
# generally intertia M = J*w_o; for a hydro turbine of 1 MW, J = 50 kg.m^2
# for 500 MW hydro, M = 50*2*pi*f (we will try to balance variation of 1500-2000 MW production from solar and wind)
M = 500*314; # 157000 kg
A = [0 1;0 0]
B = [0;1/M]
G = [0 0 0; -1/M 1/M 1/M]
C = [0 1]; # should vary between 49 to 51 Hz that is 0.2*314-1.2*314 rad/sec
# Disturbances:
N = 100;
# Load forcasted
Pl = zeros(100)
Pl[1:20] .= 1;Pl[21:40] .= 0.5;Pl[41:60] .= 1;Pl[61:80] .= 0.5;Pl[81:100] .= 1
Pl = Pl.*1500
plot(Pl,label=L"P_\ell-\mathrm{load}",lc=:blue,lw=1.5)
plot!(xlabel=L"k",
framestyle=:box,ylabel=L"P[\mathrm{MW}]",legend=:topright,xlim=(0,100),ylim=(500,1800))
# Power from PV forcasted looks like parabola
x = -100/2+1:100/2
y = @. ((-x^2)/50^2+1)*(1010-1000)+1000+rand()
Ppv = y
plot(Ppv,label=L"P_\mathrm{pv}-\mathrm{solar}",lc=:red,lw=1.5)
plot!(xlabel=L"k",
framestyle=:box,ylabel=L"P[\mathrm{MW}]",legend=:topright,xlim=(0,100),ylim=(990,1020))
# Power from wind looks like distorted sinusoidal
x = 0:0.01:2*pi
length(x)
x = range(0,length=100;stop=10*pi)
y = @. (sin(x))*(1000-500)*rand()+500
Pw = y
plot(Pw)
plot(Pw,label=L"P_\mathrm{w}-\mathrm{wind}",lc=:green,lw=1.5)
plot!(xlabel=L"k",
framestyle=:box,ylabel=L"P[\mathrm{MW}]",legend=:topright,xlim=(0,100))#,ylim=(990,1020))
# Defining reference profile-angular frequency
rN = 314*ones(N)
plot(rN,lw=1.5,ls=:solid,lc=:green,label=L"r_k",xlim=(0,100),ylim=(300,320),framestyle=:box)
# stochastic values of load solar and wind
wN = [Pl Ppv Pw]
wN[1,:] # accessing values
I tried using Optim.jl to find the optimal u_k i.e P_{m,k}. The code looks like this,
# Defining performance horizon
Np = 20
# Defining criterion function
p = 0.001
function criterion(k,xk,U,p)
R = rN[k:k+Np-1]
W = wN[k:k+Np-1,:]
J = 0.0
xi = xk
for i in 1:Np
# This is the controller model
xip1 = A*xi + B*U[i] + G*W[i,:]
xi = xip1
yi = C*xi
J = J + (yi-R[i:i])[1]^2 + p*U[i:i][1]^2
end
return J
end
;
x0 = [0.,0.]
Yc = zeros(N+1-Np)
Uc = zeros(N+1-Np)
#
Yc[1:1] = C*x0
#
xi = x0
for i in 1:(N-Np)
U0 = zeros(Np)
# assume xi is known if not estimate by Kalman filter
cost = U -> criterion(i,xi,U,p)
sol = optimize(cost,U0,NelderMead())
Uc[i] = Optim.minimizer(sol)[1]
# send control signal to real plant and collect outputs
# A, B, G, and C are different in controller and real plants
# The line below is the real system
xip1 = A*xi + B*Uc[i] + G*wN[i,:]
xi = xip1
yi = C*xi
# store real outputs
Yc[i+1:i+1] .= yi
end
plot(Yc,lw=2,lc=:blue,label=L"y_k")
plot!(rN,lw=2,ls=:solid,lc=:green,label=L"r_k")
plot!(xlabel=L"time index $k$",ylabel=L"\omega\ [rad/s]",framestyle=:box,
title="Controlled output vs. reference, p=$p",xlim=(1,80))
plot(Yc,lw=2,lc=:blue,label=L"P_{\mathrm{m},k}")
plot!(xlabel=L"time index $k$",ylabel=L"P\ [MW]",framestyle=:box,
title="Mechanical power from hydro (Controlled input), p=$p",xlim=(1,80))
The angular freq. and mechanical power looks like below from Optim.jl, however, I can not put constraint on angular frequency.
I first tried with JuMP but the optimal values it gives is always 314 rad/s for angular freq. and 0 MW for mechanical power.
The JuMP test code is like this,
using JuMP
using Ipopt
m = Model(Ipopt.Optimizer)
Np = 20;
p = 0.001
y_ini = 314*rand(Np);
u_ini = 1000*rand(Np);
@variable(m, y[i=1:Np], start=y_ini[i])
@variable(m, u[i=1:Np], start=u_ini[i])
@NLobjective(m, Min, sum((y[i]-314)^2+p*u[i]^2 for i=1:Np))
@NLconstraint(m, [i = 1:Np], 0.98*314 <= y[i] <= 1.02*314)
JuMP.optimize!(m)
JuMP.value.(u) #gives 0 for all 20 values.
JuMP.value.(y) #gives 314 for all 20 values.
My understanding of control engineering is a bit poor and I am in the learning phase and I posted this whole problem because it can be helpful for beginners like me.