SIMULATION 1
# Inspired by https://www.interviewqs.com/ddi-code-snippets/dynamic-update-plot-loop-ipython
%matplotlib inline
import time
import pylab as pl
from IPython import display

#Level Control

Fin = 6 #m3/s
Fout = 5 #m3/s
A = 3 #m2
hIC = 4 #m, initial value of the level
timePlot = 0 #s, initial time
hOld = hIC

#Integration step for Explicit Euler
intStep = 0.001;

timeArray = [timePlot]
hArray = [hIC]

pl.plot(timeArray,hArray)
pl.xlim(0,0.1)
pl.ylim(3.99,4.06)
pl.xlabel('Time (s)')
pl.ylabel('Level (m)')
display.clear_output(wait=True)
display.display(pl.gcf())
time.sleep(0.05)

for i in range(100):
    hNew = hOld + intStep*(Fin - Fout)/A
    timePlot = timePlot + intStep
    
    timeArray = pl.append(timeArray,timePlot)
    hArray = pl.append(hArray,hNew)
    
    pl.plot(timeArray,hArray)
    pl.xlim(0,0.1)
    pl.ylim(3.99,4.06)
    pl.xlabel('Time (s)')
    pl.ylabel('Level (m)')
    display.clear_output(wait=True)
    display.display(pl.gcf())
    time.sleep(0.05)
    hOld = hNew





 
SIMULATION 2 %Change tauI to 0.01 and the time at the end to 300
%matplotlib inline
import time
import pylab as pl
from IPython import display

#Level Control

Fin = 6 #m3/s
Fouts = 5 #m3/s, initial/steady-state value of Fout
A = 3 #m2
hIC = 4 #m, initial value of the level
timePlot = 0 #s, initial time
hOld = hIC
zetaIC = 0 #initial value of the controller's dynamic state
zetaOld = zetaIC
hsp = hIC #m, level set-point

Kc = -1000
tauI = 100

#Integration step for Explicit Euler
intStep = 0.0001;

timeArray = [timePlot]
hArray = [hIC]
timeArrayForu =[]
FoutArray =[]

pl.figure(1)
pl.plot(timeArray,hArray)
pl.xlim(0,0.1)
pl.ylim(3.995,4.005)
pl.xlabel('Time (s)')
pl.ylabel('Level (m)')
display.clear_output(wait=True)
display.display(pl.gcf())
time.sleep(0.001)

for i in range(200):

    Fout = Fouts + Kc*(hsp - hOld) + Kc/tauI*zetaOld
    if (Fout < 0):
        Fout = 0
        
    
    zetaNew = zetaOld + intStep*(hsp - hOld)
    
    FoutArray = pl.append(FoutArray,Fout)
    timeArrayForu = pl.append(timeArrayForu,timePlot)
   
    
    hNew = hOld + intStep*(Fin - Fout)/A
    timePlot = timePlot + intStep
    
    timeArray = pl.append(timeArray,timePlot)
    hArray = pl.append(hArray,hNew)
    
    pl.figure(1)
    pl.plot(timeArray,hArray)
    pl.xlim(0,0.1)
    pl.ylim(3.995,4.005)
    pl.xlabel('Time (s)')
    pl.ylabel('Level (m)')
    display.clear_output(wait=True)
    display.display(pl.gcf())
    time.sleep(0.001)
    hOld = hNew
    zetaOld = zetaNew
    

pl.figure(2)
pl.plot(timeArrayForu,FoutArray)
pl.xlim(0,0.1)
pl.ylim(0,7)
pl.xlabel('Time (s)')
pl.ylabel('Fout (m3/s)')
 
SIMULATION 3
%matplotlib inline
import time
import pylab as pl
from IPython import display


#xdot = x^2 + u
#u = -x^2 - x
#V = x^2

intStep = 0.001 #integration step for Explicit Euler
xIC = 2 #initial value of x
xOld = xIC #First value of xOld
curTime = 0

uArray = []
timeArrayForu = []
timeArray = [curTime]
VArray =[xOld**2]
xArray =[xOld]

pl.figure(1)
pl.plot(timeArray,VArray)
pl.xlim(0,0.1)
pl.ylim(0,5)
pl.xlabel('Time (s)')
pl.ylabel('V')
display.clear_output(wait=True)
display.display(pl.gcf())
time.sleep(0.001)


for i in range(100):
    u = -xOld**2 - xOld
    xNew = xOld + intStep*(xOld**2 + u)
    
    uArray = pl.append(uArray,u)
    timeArrayForu = pl.append(timeArrayForu,curTime)

    curTime = curTime + intStep
    xOld = xNew
    
    timeArray = pl.append(timeArray,curTime)
    VArray = pl.append(VArray,xNew**2)
    xArray = pl.append(xArray,xNew)
    
    pl.figure(1)
    pl.plot(timeArray,VArray)
    pl.xlim(0,0.1)
    pl.ylim(2.5,4.5)
    pl.xlabel('Time (s)')
    pl.ylabel('V')
    display.clear_output(wait=True)
    display.display(pl.gcf())
    time.sleep(0.001)

pl.figure(2)
pl.plot(timeArrayForu,uArray)
pl.xlim(0,0.1)
pl.ylim(-7,0)
pl.xlabel('Time (s)')
pl.ylabel('u')

pl.figure(3)
pl.plot(timeArray,xArray)
pl.xlim(0,0.1)
pl.ylim(1,2.5)
pl.xlabel('Time (s)')
pl.ylabel('x')



 
SIMULATION 4   https://colab.research.google.com/github/ndcbe/CBE60499/blob/master/docs/01.00-Pyomo-Introduction.ipynb#scrollTo=hvBvTJBMhbDT
https://ndcbe.github.io/CBE60499/01.00-Pyomo-Introduction.html
https://github.com/ndcbe/CBE60499/blob/main/notebooks/helper.py

import sys
import os

os.system("pip install -q idaes_pse")
os.system('idaes –version')
os.system("idaes get-extensions")
os.system("ln -s /root/.idaes/bin/ipopt ipopt")

os.system('wget -N -q "https://ampl.com/dl/open/ipopt/ipopt-linux64.zip"')
os.system('!unzip -o -q ipopt-linux64')

import pyomo.environ as py
pyModel = py.ConcreteModel()
pyModel.x1 = py.Var(initialize=1, bounds=(-5,5))
pyModel.x2 = py.Var(initialize=1, bounds=(0,10))
pyModel.x3 = py.Var(initialize=1, bounds=(-1,2))

pyModel.objectiveFtn = py.Objective(expr=2*pyModel.x1+3*pyModel.x2**2-5*pyModel.x3, sense=py.minimize)
pyModel.con1 = py.Constraint(expr= 3*pyModel.x1 + 2*pyModel.x2 == 15)
pyModel.con2 = py.Constraint(expr = 5*pyModel.x1 + pyModel.x3 == 10)

pyModel.pprint()


opt1 = py.SolverFactory('ipopt')
status1 = opt1.solve(pyModel, tee=True)


print("x1 = ",py.value(pyModel.x1))
print("x2 = ",py.value(pyModel.x2))
print("x3 = ",py.value(pyModel.x3))
print("\n")


import numpy as np
import matplotlib.pyplot as pl

x3_try = np.linspace(-1,0.1,2)
x3_tryAllowed = []

objVector = []

for x3Val in x3_try:
  A = np.array([[3, 2], [5, 0]])
  b = np.array([15, 10-x3Val])
  xsoln = np.linalg.solve(A,b)
  x1Val = xsoln[0]
  x2Val = xsoln[1]
  if x1Val >= -5 and x1Val <= 5:
    if x2Val >= 0 and x2Val <= 10:
      objVal = 2*x1Val+3*x2Val**2-5*x3Val
      objVector=np.append(objVector,objVal)
      x3_tryAllowed = np.append(x3_tryAllowed,x3Val)

pl.plot(x3_tryAllowed,objVector,color='black',linewidth=2,label="Test Cases")
pl.xlabel("$x_3$",fontsize = 20)
pl.ylabel("Objective Function Value",fontsize=20)

x3_optimal = py.value(pyModel.x3)
obj_optimal = py.value(pyModel.objectiveFtn)
pl.plot(x3_optimal,obj_optimal,marker='*',color='blue',markersize=20,label="Optimal Value")
pl.legend()
pl.grid(True)
pl.show()

 
SIMULATION 5


import sys
import os

os.system("pip install -q idaes_pse")
os.system('idaes –version')
os.system("idaes get-extensions")
os.system("ln -s /root/.idaes/bin/ipopt ipopt")

os.system('wget -N -q "https://ampl.com/dl/open/ipopt/ipopt-linux64.zip"')
os.system('!unzip -o -q ipopt-linux64')

import pyomo.environ as py
import math
import numpy as np

hint = 0.01 #integration step for explicit Euler
Delta = 0.02 #sampling period length
N = 2 #Number of sampling periods in prediction horizon
numSampTotal = 2 #Number of sampling periods to be simulated in the simulation

xinitial = 2 #Value of x(t0) at the initial time
xOld = xinitial
timePlot = 0;
xVector =[xOld]
uVector = []
timeVector = [timePlot]
timeForuVector =[]

for i in range(numSampTotal):
  #Call decision variables x(t1), x(t2), x(t3), and x(t4) as x1, x2, x3, and x4, respectively

  pyModel = py.ConcreteModel()
  pyModel.x1 = py.Var(initialize=1, bounds=(-3,3))
  pyModel.x2 = py.Var(initialize=1, bounds=(-3,3))
  pyModel.x3 = py.Var(initialize=1, bounds=(-3,3))
  pyModel.x4 = py.Var(initialize=1, bounds=(-3,3))

  #Call decision variables u(t0) and u(t1) as x5 and x6, respectively
  pyModel.x5 = py.Var(initialize=1, bounds=(-1,1))
  pyModel.x6 = py.Var(initialize=1, bounds=(-1,1))
  
  pyModel.objectiveFtn = py.Objective(expr=hint*(7*pyModel.x1 - 3*pyModel.x1**2)+hint*(7*pyModel.x2 - 3*pyModel.x2**2) + hint*(7*pyModel.x3 - 3*pyModel.x3**2)+ hint*(7*pyModel.x4 - 3*pyModel.x4**2), sense=py.minimize)
  pyModel.con1 = py.Constraint(expr= -pyModel.x1 + xOld + hint*(2*xOld**2 + xOld + pyModel.x5)  == 0)
  pyModel.con2 = py.Constraint(expr= -pyModel.x2 + pyModel.x1 + hint*(2*(pyModel.x1)**2 + pyModel.x1 + pyModel.x5)  == 0)
  pyModel.con3 = py.Constraint(expr= -pyModel.x3 + pyModel.x2 + hint*(2*(pyModel.x2)**2 + pyModel.x2 + pyModel.x6)  == 0)
  pyModel.con4 = py.Constraint(expr= -pyModel.x4 + pyModel.x3 + hint*(2*(pyModel.x3)**2 + pyModel.x3 + pyModel.x6)  == 0)

  opt1 = py.SolverFactory('ipopt')
  status1 = opt1.solve(pyModel, tee=True)

  
  #Store the first input in the prediction horizon
  u1 = py.value(pyModel.x5)

  uVector = np.append(uVector,u1)
  uVector = np.append(uVector,u1)
  timeForuVector = np.append(timeForuVector,timePlot)
  timeForuVector = np.append(timeForuVector,timePlot+Delta)

  #Apply the first input for a sampling period (we designed this problem to have 2 integration steps in a sampling period)
  for j in range(2):
    xNew = xOld + hint*(2*xOld**2 + xOld + u1)
    xOld = xNew
    timePlot = timePlot + hint
    xVector = np.append(xVector,xNew)
    timeVector = np.append(timeVector,timePlot)


import matplotlib.pyplot as pl

pl.figure(1)
pl.plot(timeVector,xVector,color='black',linewidth=2)
pl.xlabel("Time",fontsize = 20)
pl.ylabel("$x$",fontsize=20)
pl.show()

pl.figure(2)
pl.plot(timeForuVector,uVector,color='black',linewidth=2)
pl.xlabel("Time",fontsize = 20)
pl.ylabel("$u$",fontsize=20)
pl.show()

print(u1)
print(uVector)