Dynamic Optimization

ax.plot(t, p,'b:', label='Production 1 ($g_1$)',\
        linewidth=2)

ax.plot(t,recovery,'b:',label='Recovered ($e_{\text{out},1}$)',\
        linewidth=2)

       label='Production 2 ($g_2$)',linewidth=2)

       label='Recovered ($e_{\text{out},2}$)',linewidth=2)

* Gates, N.S., Hill, D.C., Billings, B.W., Powell, K.M., Hedengren, J.D., Benchmarks for Grid Energy Management with Python Gekko, 60th IEEE Conference on Decision and Control (CDC), Tutorial Session: Open Source Software for Control Applications, Austin, TX, Dec 13-15, 2021. [[Attach:Gates_2021.pdf|Preprint]]

e  = m.CV(0)
e.STATUS=1; e.SPHI=e.SPLO=0
e.WSPHI=1000; e.WSPLO=1

m.Equations([g.dt()==r, e==d-g])

plt.plot(t,g,'r-',label='Production')
plt.plot(t,d,'b:',label='Demand')

e1  = m.CV(0); e1.STATUS=1; e1.SPHI=e1.SPLO=0; e1.WSPHI=1000; e1.WSPLO=1
e2  = m.CV(0); e2.STATUS=1; e2.SPHI=e2.SPLO=0; e2.WSPHI=1000; e2.WSPLO=1

m.Equations([g1.dt()==r, e1==d1-g1, e2==d2-g2])

plt.plot(t,d1,'r-',label='Prod 1')
plt.plot(t,g1,'b:',label='Demand 1')

plt.plot(t,d2,'r-',label='Prod 2')
plt.plot(t,g2,'b:',label='Demand 2')

<iframe width="560" height="315" src="https://www.youtube.com/embed/nsBeSUBzurc" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>

<iframe width="560" height="315" src="https://www.youtube.com/embed/7OhPDdioke8" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>

%width=550px%Attach:grid_benchmarks.png

Source code examples from the publication:

g1  = m.Var(d1[0]); dh = m.Intermediate(g1*2)

e  = m.CV(0); e.STATUS=1; e.SPHI=e.SPLO=0; e.WSPHI=1000; e.WSPLO=1
h  = m.CV(0); h.STATUS=1; h.SPHI=h.SPLO=0; h.WSPHI=1000; h.WSPLO=1

m.Equations([g1.dt()==r, e==d1-g1, h==d2-dh])

plt.plot(t,dh,'b:',label='Demand 2')

The Gekko Optimization Suite is a machine learning and optimization package in Python for mixed-integer and differential algebraic equations and is capable of solving complex grid design and control problems. A series of non-dimensional benchmark cases are proposed for grid energy production. These include:

fe = m.Param(np.cos(2*np.pi*t)+3)
fh = m.Param(1.5*np.sin(2*np.pi*t)+7)
fz = m.Param(np.clip(-0.2*np.sin(2*np.pi*t),0,None))
de = m.Var(fe[0]); dh=m.Intermediate(de*2)
dz = m.Var(0,lb=0); dze=m.Var(0); dzh=m.Var(0)
te = m.Intermediate(fe+dze); th = m.Intermediate(fh+dzh)

der = m.MV(0,lb=-1,ub=1); der.STATUS=1; der.DCOST=0.0
dzr = m.MV(0,lb=-1,ub=1); dzr.STATUS=1; dzr.DCOST=0.0

m.Equations([de.dt()==der, e==te-de, h==th-dh])
m.Equations([dz.dt()==dzr, z==fz-dz, dze==dz*2, dzh==dz*3])
m.Maximize(dz)

plt.plot(t,te,'r-',label='Demand 1')
plt.plot(t,de,'b:',label='Prod 1')

plt.plot(t,th,'r-',label='Demand 2')
plt.plot(t,dh,'b:',label='Prod 2')

plt.plot(t,fz,'r-',label='Demand 3')
plt.plot(t,dz,'b:',label='Prod 3')

plt.plot(t,der,'k:',label='Ramp Rate 1')
plt.plot(t,dzr,'k--',label='Ramp Rate 3')

fe  = m.Param(np.cos(2*np.pi*t)+3)
fh  = m.Param(1.5*np.sin(2*np.pi*t)+7)
de  = m.Var(fe[0]); dh = m.Intermediate(de*2)

der = m.MV(0,lb=-1,ub=1); der.STATUS=1

m.Equations([de.dt()==der, e==fe-de, h==fh-dh])

plt.plot(t,fe,'r-',label='Prod 1')
plt.plot(t,de,'b:',label='Demand 1')

plt.plot(t,fh,'r-',label='Prod 2')

plt.plot(t,der,'k--',label='Ramp Rate')

fe  = m.Param(np.cos(2*np.pi*t)+3)
de  = m.Var(fe[0])

der = m.MV(0,lb=-1,ub=1); der.STATUS=1

m.Equations([de.dt()==der, e==fe-de])

plt.plot(t,fe,'r-',label='Production')
plt.plot(t,de,'b:',label='Demand')
plt.plot(t,der,'k--',label='Ramp Rate')

Source code examples from the publication (under review):

* Gates, N.S., Hill, D.C., Billings, B.W., Powell, K.M., Hedengren, J.D., Benchmarks for Grid Energy Management with Python Gekko, 60th IEEE Conference on Decision and Control (CDC), Tutorial Session: Open Source Software for Control Applications, Austin, TX, Dec 13-15, 2021.

{$ $$\begin{align} \min_g \ \ & g \\ \mathrm{s.t.} & \frac{de}{dt} = e_{in} - e_{out} \cdot \eta \\ & e_{in} = g - d + s_{in} \\ & e_{out} = d - g + s_{out} \\ & g - d = s_{out} - s_{in} \\ & s_{out}, \, s_{in} \geq 0, \enspace e_{out} \cdot e_{in} \leq 0 \\ & g + e_{out}/\eta - e_{in} \geq d \\ & e \geq 0, \enspace \eta = 0.7 \\ & d = 10 - 2 \sin(2 \, \pi \, t) \\ & e(0) = e(1) = 0 \end{align}$$ $} 

{$ $$\begin{align} \min_{r} \ \ & J = \sum_{i=1}^n \int_{t=0}^{1} \, [ 1000 \, \max(0, \, \Psi) \vphantom{\int} + \, \max(0, \, -\Psi) \vphantom{\int} ]dt \\ & \text{where} \ \Psi = d - g - R + e_{out}/\eta - e_{in} \\ \mathrm{s.t.}\ \ & \frac{de}{dt} = e_{in} - e_{out} \cdot \eta \\ & e_{in} = g + R - d + s_{in} \\ & e_{out} = d - g - R + s_{out} \\ & g + R - d = s_{out} - s_{in} \\ & s_{out}, \, s_{in} \geq 0, \enspace e_{out} \cdot e_{in} \leq 0\\ & g + R + e_{out}/\eta - e_{in} \geq d \\ & e \geq 0, \enspace \eta = 0.85 \\ & d = 7 - 2 \sin(2 \, \pi \, t) \\ & R = \begin{cases}  3 + 3 \cos(4 \, \pi \, t) & \frac{1}{4} \leq t \leq \frac{3}{4} \\ 0 & \mbox{otherwise} \\ \end{cases} \\ & \frac{dg}{dt} = r, \enspace -1\leq r \leq 1 \\ & e(0) = e(1) = 0 \end{align}$$ $}

{$ \begin{align} \min_{r_1} \ \ & g_1 \\ \mathrm{s.t.}\ \ & \frac{de_i}{dt} = e_{in,i} - e_{out,i}\cdot \eta_i \\ & e_{in,i} = g_i + R_i - d_i + s_{in,i} \\ & e_{out,i} = d_i - g_i - R_i + s_{out,i} \\ & g_i + R_i - d_i = s_{out,i} - s_{in,i} \\ & s_{out,i}, \, s_{in,i} \geq 0, \enspace e_{out,i} \cdot e_{in,i} \leq 0 \\ & g_i + R_i + e_{out,i}/\eta_i - e_{in,i} \geq d_i \\ & e_i \geq 0, \enspace \eta_1 = 0.7, \enspace \eta_2 = 0.8 \\ & d_1 = 10 - 2 \sin(2 \, \pi \, t) \\ & d_2 = 15 + 1.5 \cos(2 \, \pi \, t) \\ & R_1 = $$\begin{cases} 2 + 2 \cos(4 \, \pi \, t) & \frac{1}{4} \leq t \leq \frac{3}{4} \\ 0 & \mbox{otherwise} \\ \end{cases}$$ \\ & \frac{dg_1}{dt} = r_1, \enspace -3\leq r_1 \leq 3 \\ & R_2 = 0, \enspace g_2 = 1.5 \cdot g_1 \\ & e_1(0) = 0, \enspace e_2(0) = e_2(1) = 0.5 \end{align} $}

ax.plot(t,stored,'g--',label='Stored ($e_{in,1}$)')
ax.plot(t,recovery,'b:',label='Recovered ($e_{out,1}$)',\

ax.plot(t,stored_h,'g--',label='Stored ($e_{in,2}$)')

       label='Recovered ($e_{out,2}$)',linewidth=2)

|| {`\Phi`} || objective function ||

{$ \begin{align}    \min_{r} \    & \Phi = \int_{t=0}^{1} \left [ 1000 \max(0, d - g) + \max(0, g - d) \right ] \\ \mathrm{s.t.}\  & \frac{dg}{dt} = r \\ & d=\cos\left({2\pi t}\right) + 3 \label{eqn:benchmark1c}\\ & -1 \leq r \leq 1 \label{eqn:benchmark1d}\\ & g(0)=d(0)=4 \end{align} $}

{$ \begin{align} \min_{r} \ \ & \Phi = \sum_{i=1}^n \int_{t=0}^{1} \, \left [ 1000 \; \max(0, \, d_i - g_i) \vphantom{\int} + \max(0, \, g_i - d_i) \vphantom{\int} \right ] \\ \mathrm{s.t.}\ \ & \frac{dg_1}{dt} = r, \enspace g_2 = 2 g_1 \label{eqn:benchmark2b}\\ & d_1=\cos\left({2\pi t}\right) + 3 \\ & d_2=1.5\sin\left({2\pi t}\right) + 7 \\ & -1 \leq r \leq 1 \\ & g_1(0)=d_1(0)=4 \\ & g_2(0)=8, \; d_2(0)=7 \end{align} $}

{$ \begin{align} \min_{r_1, r_3} \ \ & \Phi = \sum_{i=1}^n \int_{t=0}^{1} \left [ \vphantom{\int} 1000 \; \max(0, \, d_i - g_i) + \max(0, \, g_i - d_i) \vphantom{\int} \right ] - 0.1 \int_{t=0}^{1} g_3 \\  \mathrm{s.t.}\ \ & \frac{dg_1}{dt} = r_1, \; \frac{dg_3}{dt} = r_3 \\ & g_2 = 2 g_1 \\ & d_1=\cos\left({2\pi t}\right) + 3 + 2 g_3 \\  & d_2=1.5\sin\left({2\pi t}\right) + 7 + 3 g_3 \\  & d_3=\max\left[0,-0.2\sin\left(2\pi t\right) \right] \\ & -1 \leq r_1 \leq 1, \enspace -1 \leq r_3 \leq 1 \\ & g_1(0)=d_1(0)=4 \\ & g_2(0)=8, \; d_2(0)=7 \\  & g_3(0)=d_3(0)=0 \end{align} $}

{$ \begin{align} \min_{r} \ \ & \Phi = \sum_{i=1}^n \int_{t=0}^{1} \, [ 1000 \, \max(0, \, \Psi) \vphantom{\int} + \, \max(0, \, -\Psi) \vphantom{\int} ] \\ & where \ \Psi = d - g - R + e_{out}/\eta - e_{in} \\ \mathrm{s.t.}\ \ & \frac{de}{dt} = e_{in} - e_{out} \cdot \eta \\ & e_{in} = g + R - d + s_{in} \\ & e_{out} = d - g - R + s_{out} \\ & g + R - d = s_{out} - s_{in} \\ & s_{out}, \, s_{in} \geq 0, \enspace e_{out} \cdot e_{in} \leq 0\\ & g + R + e_{out}/\eta - e_{in} \geq d \\ & e \geq 0, \enspace \eta = 0.85 \\ & d = 7 - 2 \sin(2 \, \pi \, t) \\ & R = $$\begin{cases}  3 + 3 \cos(4 \, \pi \, t) & \frac{1}{4} \leq t \leq \frac{3}{4} \\ 0 & \mbox{otherwise} \\ \end{cases}$$ \\ & \frac{dg}{dt} = r, \enspace -1\leq r \leq 1 \\ & e(0) = e(1) = 0 \end{align} $}

             err == d - g - r + recover/eta - store,                                

             s.dt() == store - recover/eta,                

p = m.FV(); p.STATUS = 1 # production

vy = m.Var(lb=0)        # store slack variable

vx = m.Var(lb=0)        # recover slack variable

m.Equations([p + recover/eta - store >= d,
            p - d == vx - vy,
            store == p - d + vy,
            recover == d - p + vx,

m.Minimize(p)

plt.plot(m.time,p,'b:',label='Prod')

Source code examples from the publication (to appear): Benchmarks for Grid Energy Management with Python Gekko

der = m.MV(0,lb=-1,ub=1); der.STATUS=1; der.DCOST=0.01
dzr = m.MV(0,lb=-1,ub=1); dzr.STATUS=1; dzr.DCOST=0.01

plt.subplot(4,1,1); plt.plot(t,te,'r-'); plt.plot(t,de,'b:'); plt.grid()
plt.subplot(4,1,2); plt.plot(t,th,'r-'); plt.plot(t,dh,'b:'); plt.grid()
plt.subplot(4,1,3); plt.plot(t,fz,'r-'); plt.plot(t,dz,'b:'); plt.grid()
plt.subplot(4,1,4); plt.plot(t,der,'k:'); plt.plot(t,dzr,'k--'); plt.grid()
plt.xlabel('Time'); plt.show()

!!!! Benchmark IV: Constant Production with Energy Storage

!!!! Benchmark V: Load Following with Energy Storage

!!!! Benchmark VI: Cogeneration with Dual Energy Storage

{$ \begin{align} \min_g \ \ & g \\ \mathrm{s.t.} & \frac{de}{dt} = e_{in} - e_{out} \cdot \eta \\ & e_{in} = g - d + s_{in} \\ & e_{out} = d - g + s_{out} \\ & g - d = s_{out} - s_{in} \\ & s_{out}, \, s_{in} \geq 0, \enspace e_{out} \times e_{in} \leq 0 \\ & g + e_{out}/\eta - e_{in} \geq d \\ & e \geq 0, \enspace \eta = 0.7 \\ & d = 10 - 2 \sin(2 \, \pi \, t) \\ & e(0) = e(1) = 0 \end{align} $} 

{$ \begin{align} \min_{r} \ \ & \Phi = \sum_{i=1}^n \int_{t=0}^{1} \, [ 1000 \, \max(0, \, \Psi) \vphantom{\int} + \, \max(0, \, -\Psi) \vphantom{\int} ] \\ & where \ \Psi = d - g - R + e_{out}/\eta - e_{in} \\ \mathrm{s.t.}\ \ & \frac{de}{dt} = e_{in} - e_{out} \cdot \eta \\ & e_{in} = g + R - d + s_{in} \\ & e_{out} = d - g - R + s_{out} \\ & g + R - d = s_{out} - s_{in} \\ & s_{out}, \, s_{in} \geq 0, \enspace e_{out} \times e_{in} \leq 0\\ & g + R + e_{out}/\eta - e_{in} \geq d \\ & e \geq 0, \enspace \eta = 0.85 \\ & d = 7 - 2 \sin(2 \, \pi \, t) \\ & R = $$\begin{cases}  3 + 3 \cos(4 \, \pi \, t) & \frac{1}{4} \leq t \leq \frac{3}{4} \\ 0 & \mbox{otherwise} \\ \end{cases}$$ \\ & \frac{dg}{dt} = r, \enspace -1\leq r \leq 1 \\ & e(0) = e(1) = 0 \end{align} $}

{$ \begin{align} \min_g \ \ & g \\ \mathrm{s.t.}\ \ & \frac{de}{dt} = e_{in} - e_{out} \cdot \eta \\ & e_{in} = g - d + s_{in} \\ & e_{out} = d - g + s_{out} \\ & g - d = s_{out} - s_{in} \\ & s_{out}, \, s_{in} \geq 0, \enspace e_{out} \times e_{in} \leq 0 \\ % & e_{out} \times e_{in} \leq 0 \\ & g + e_{out}/\eta - e_{in} \geq d \\ & e \geq 0, \enspace \eta = 0.7 \\ & d = 10 - 2 \sin(2 \, \pi \, t) \\ & e(0) = e(1) = 0 \end{align} $}

       \begin{align}
          \min_g \ \     & g \label{eqn:benchmark4a}\\
           \mathrm{s.t.}\ \  & \frac{de}{dt} = e_{in} - e_{out} \cdot \eta \label{eqn:benchmark4b} \\
                         & e_{in} = g - d + s_{in} \label{eqn:benchmark4c} \\
                          & e_{out} = d - g + s_{out} \label{eqn:benchmark4d}\\
                            & g - d = s_{out} - s_{in} \label{eqn:benchmark4e}\\
                            & s_{out}, \, s_{in} \geq 0, \enspace e_{out} \times e_{in} \leq 0 \label{eqn:benchmark4f}\\
                            % & e_{out} \times e_{in} \leq 0 \label{eqn:benchmark4g} \\
                            & g + e_{out}/\eta - e_{in} \geq d \label{eqn:benchmark4h} \\
                            & e \geq 0, \enspace \eta = 0.7 \label{eqn:benchmark4i} \\
                            & d = 10 - 2 \sin(2 \, \pi \, t) \label{eqn:benchmark4j} \\
                            & e(0) = e(1) = 0 \label{eqn:benchmark4k}
        \end{align}

The third benchmark problem enhances the previous problems further still, creating a tri-generation system ($n=3$) with two producers, three products, and three demand profiles. The primary producer is the same as the prior benchmark and is ramp-rate constrained to produce the two primary products (e.g., electricity and heat). An additional producer (e.g., a solid oxide electrolysis cell) uses these first two products to make a third product (e.g., hydrogen), thereby utilizing any excess system capacity and maximizing its production while avoiding supply shortages for products one and two. As before, the optimizer has perfect foresight.

{$ \begin{align} \min_{r} \ \ & \Phi = \sum_{i=1}^n \int_{t=0}^{1} \, \left [ 1000 \; \max(0, \, d_i - g_i) \vphantom{\int} \right. + \max(0, \, g_i - d_i) \vphantom{\int} \right ] \\ \mathrm{s.t.}\ \ & \frac{dg_1}{dt} = r, \enspace g_2 = 2 g_1 \label{eqn:benchmark2b}\\ & d_1=\cos\left({2\pi t}\right) + 3 \\ & d_2=1.5\sin\left({2\pi t}\right) + 7 \\ & -1 \leq r \leq 1 \\ & g_1(0)=d_1(0)=4 \\ & g_2(0)=8, \; d_2(0)=7 \end{align} $}

{$ \begin{align} \min_{r_1, r_3} \ \ & \Phi = \sum_{i=1}^n \int_{t=0}^{1} \left [ \vphantom{\int} 1000 \; \max(0, \, d_i - g_i) + \max(0, \, g_i - d_i) \vphantom{\int} \right ] - 0.1 \int_{t=0}^{1} g_3 \\  \mathrm{s.t.}\ \ & \frac{dg_1}{dt} = r_1, \; \frac{dg_3}{dt} = r_3 \\ & g_2 = 2 g_1 \\ & d_1=\cos\left({2\pi t}\right) + 3 + 2 g_3 \\  & d_2=1.5\sin\left({2\pi t}\right) + 7 + 3 g_3 \\  & d_3=\max\left[0,-0.2\sin\left(2\pi t\right) \right] \\  %& dt_1 = d_1 + 2 \, g_3 \\ %& dt_2 = d_1 + 3 \, g_3 \\ & -1 \leq r_1 \leq 1, \enspace -1 \leq r_3 \leq 1 \\ & g_1(0)=d_1(0)=4 \\ & g_2(0)=8, \; d_2(0)=7 \\  & g_3(0)=d_3(0)=0 \end{align} $}

        \begin{align}
           \min_{r_1, r_3} \ \   & \Phi = \sum_{i=1}^n \int_{t=0}^{1} \left [ \vphantom{\int} 1000 \; \max(0, \,  d_i - g_i) \; \label{eqn:benchmark3a} \right. \\
           & \ \ \ \ \ \ \ \left. + \max(0, \, g_i - d_i) \vphantom{\int} \right ] - 0.1 \int_{t=0}^{1} g_3 \notag \\
           \mathrm{s.t.}\ \  & \frac{dg_1}{dt} = r_1, \; \frac{dg_3}{dt} = r_3 \label{eqn:benchmark3b} \\
                          & g_2 = 2 g_1 \label{eqn:benchmark3c} \\
                            & d_1=\cos\left({2\pi t}\right) + 3 + 2 g_3 \label{eqn:benchmark3d}\\
                            & d_2=1.5\sin\left({2\pi t}\right) + 7 + 3 g_3 \label{eqn:benchmark3e}\\
                            & d_3=\max\left[0,-0.2\sin\left(2\pi t\right) \right] \label{eqn:benchmark3f}\\
                            %& dt_1 = d_1 + 2 \, g_3 \label{eqn:benchmark3h} \\
                            %& dt_2 = d_1 + 3 \, g_3 \label{eqn:benchmark3i} \\
                            & -1 \leq r_1 \leq 1, \enspace -1 \leq r_3 \leq 1 \label{eqn:benchmark3g} \\
                            & g_1(0)=d_1(0)=4 \label{eqn:benchmark3h}  \\
                            & g_2(0)=8, \; d_2(0)=7 \label{eqn:benchmark3i} \\
                            & g_3(0)=d_3(0)=0 \label{eqn:benchmark3j}
        \end{align}

{$ \begin{align} \min_{r} \ \ & \Phi = \sum_{i=1}^n \int_{t=0}^{1} \, \left [ 1000 \; \max(0, \, d_i - g_i) \vphantom{\int} \right. + \left. \max(0, \, g_i - d_i) \vphantom{\int} \right ] \\ \mathrm{s.t.}\ \ & \frac{dg_1}{dt} = r, \enspace g_2 = 2 g_1 \label{eqn:benchmark2b}\\ & d_1=\cos\left({2\pi t}\right) + 3 \\ & d_2=1.5\sin\left({2\pi t}\right) + 7 \\ & -1 \leq r \leq 1 \\ & g_1(0)=d_1(0)=4 \\ & g_2(0)=8, \; d_2(0)=7 \end{align} $}

{$ \begin{align} \min_{r} \ \ & \Phi = \sum_{i=1}^n \int_{t=0}^{1} \, \left [ 1000 \; \max(0, \, d_i - g_i) \vphantom{\int} \right. \label{eqn:benchmark2a} \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ + \; \left. \max(0, \, g_i - d_i) \vphantom{\int} \right ] \\ \mathrm{s.t.}\ \ & \frac{dg_1}{dt} = r, \enspace g_2 = 2 g_1 \label{eqn:benchmark2b}\\ & d_1=\cos\left({2\pi t}\right) + 3 & d_2=1.5\sin\left({2\pi t}\right) + 7 \\ & -1 \leq r \leq 1 \\ & g_1(0)=d_1(0)=4 \\ & g_2(0)=8, \; d_2(0)=7 \end{align} $}

||border=0|
| '''Symbol''' | '''Description''' |
| {`\Phi`} | objective function |
| ''d'' | demand |
| ''g'' | generation |
| ''r'' | ramp rate |
| ''e'' | storage inventory |
| ''e'_in_''', ''e'_out_''' | energy stored, recovered |
| ''R'' | renewable source |
| ''s'_in_''', ''s'_out_''' | slack variables for storage switching |
| {`\eta`} | storage efficiency |
| ''n'' | number of generating units |
| ''i'' | subscript indicates product ''i'' |

(:keywords nonlinear control, energy, optimal control, dynamic optimization, engineering optimization, MATLAB, Python, GEKKO, differential, algebraic, modeling language, university course:)

(:title Optimal Control Benchmark Problems III:)