# Stochastic MPC using JuMP - Explained in detail for future example in JuMP

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;
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!(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)
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.

1 Like

Is there a particular question? The optimal solution for u is always 0 because you have no constraints linking y and u in your JuMP model.

I could not formulate JuMP optimization model from the cost function given above which include stochasticity of wind and solar and multi-state. How do I do that?