text: # GP NMPC problem setup
  import numpy as np
  from casadi import *
  
  
  class Bio_reactor_1:
      def __init__(self):
          self.real_parameters =  [0.0923*0.62, 178.85]
      def specifications(self):
          ''' Specify Problem parameters '''
          tf              = 240.      # final time
          nk              = 12        # sampling points
          dt              = tf/nk
          x0              = np.array([1.,150.,0.])
          Lsolver         = 'mumps'  #'ma97'  # Linear solver
          c_code          = False    # c_code
          shrinking_horizon = True
  
          return dt, x0, Lsolver, c_code, shrinking_horizon
  
      def DAE_system(self, uncertain_parameters=False):
          # Define vectors with names of states
          states     = ['x','n','q']
          nd         = len(states)
          xd         = SX.sym('xd',nd)
          for i in range(nd):
              globals()[states[i]] = xd[i]
  
          # Define vectors with names of algebraic variables
          algebraics = []
          na         = len(algebraics)
          xa         = SX.sym('xa',na)
          for i in range(na):
              globals()[algebraics[i]] = xa[i]
  
          set_points = []
          n_ref      = len(set_points)
          x_ref      = SX.sym('x_ref',n_ref)
          for i in range(n_ref):
              globals()[set_points[i]] = x_ref[i]
  
  
  
          # Define vectors with banes of input variables
          inputs     = ['L','Fn']
          nu         = len(inputs)
          u          = SX.sym("u",nu)
          for i in range(nu):
              globals()[inputs[i]] = u[i]
  
          # Define model parameter names and values
          if uncertain_parameters:
              modpar = ['u_m', 'k_s']
              nmp = len(modpar)
              theta = SX.sym("theta", nmp)
              for i in range(nmp):
                  globals()[modpar[i]] = theta[i]
          else:
              modpar = ['u_m', 'k_s']
              nmp = len(modpar)
              theta = SX.sym("theta", nmp)
              for i in range(nmp):
                  globals()[modpar[i]] = theta[i]
  
          k_i, K_N, u_d, Y_nx, k_m, k_sq,\
          k_iq, k_d, K_Np = 447.12, 393.10, 0.001, 504.49,\
                            2.544*0.62*1e-4, 23.51, 800.0, 0.281, 16.89
          # Additive measurement noise
      #    Sigma_v  = [400.,1e5,1e-2]*diag(np.ones(nd))*1e-6
  
          # Additive disturbance noise
      #    Sigma_w  = [400.,1e5,1e-2]*diag(np.ones(nd))*1e-6
  
          # Initial additive disturbance noise
      #    Sigma_w0 = [1.,150.**2,0.]*diag(np.ones(nd))*1e-3
          # Declare ODE equations (use notation as defined above)
  
          dx   = u_m * L/(L+k_s+L**2./k_i) * x * n/(n+K_N) - u_d*x
          dn   = - Y_nx*u_m* L/(L+k_s+L**2./k_i) * x * n/(n+K_N)+ Fn
          dq   = (k_m * L/(L+k_sq+L**2./k_iq) * x - k_d * q/(n+K_Np))# * (sign(499.9999 - n)+1)/2 * (sign(x - 10.0001)+1)/2
  
                # * (tanh(1000*(500. - n))+1)/2 * (tanh(1000*(x - 10.000))+1)/2# * (sign(499.9999 - n)+1)/2 * (sign(x - 10.0001)+1)/2
  
          ODEeq =  [dx, dn, dq]
  
          # Declare algebraic equations
          Aeq = []
  
          # Define objective to be minimized
          t     = SX.sym('t')
          Obj_M = 0.#Function('mayer', [xd, x_ref, u], [-q])  # Mayer term
          Obj_L = 0.#Function('lagrange', [xd, x_ref, u], [0.])  # Lagrange term
          Obj_D = 0.#Function('Discrete', [xd, x_ref, u], [0.])  # Lagrange term
  
          R           = np.diagflat([3.125e-8, 3.125e-006])                         # Weighting of control penality
  
          # Define control bounds
          u_min = np.array([120., 0.])
          u_max = np.array([400., 40.])
          x_min = np.array([0.]*nd)
          x_max = np.array([np.inf]*nd)
          # Define constraint functions g(x) <= 0
          gequation = vertcat(n - 800.)#n - 800., q - 0.011 * x)
          ng = SX.size(gequation)[0]
          gfcn = Function('gfcn', [xd, xa, u], [gequation])
  
          return xd, xa, u, theta, ODEeq, Aeq, u_min, u_max, x_min, x_max, states,\
                 algebraics, inputs, nd, na, nu, n_ref, nmp, self.real_parameters, ng, gfcn, Obj_M, Obj_L, Obj_D, R
  
  
  
      def integrator_model(self):
          """
          This function constructs the integrator to be suitable with casadi environment, for the equations of the model
          and the objective function with variable time step.
           inputs: NaN
           outputs: F: Function([x, u, dt]--> [xf, obj])
          """
          dt, x0, Lsolver, c_code, shrinking_horizon = self.specifications()
  
          xd, xa, u, uncertainty, ODEeq, Aeq, u_min, u_max, x_min, x_max, states, \
          algebraics, inputs, nd, na, nu, n_ref, nmp, modparval, ng, gfcn, Obj_M, Obj_L, Obj_D, R = self.DAE_system()
          dae = {'x': vertcat(xd), 'z': vertcat(xa), 'p': vertcat(u, uncertainty),
                 'ode': vertcat(*ODEeq), 'alg': vertcat(*Aeq)}
  
          opts = {'tf': dt}  # interval length
          F = integrator('F', 'idas', dae, opts)
  
          return F
  
  
  class Bio_reactor_2:
      def __init__(self):
          self.real_parameters =  [0.0923*0.62, 178.85]
      def specifications(self):
          ''' Specify Problem parameters '''
          tf              = 240.      # final time
          nk              = 12        # sampling points
          dt              = tf/nk
          x0              = np.array([1.,150.,0.])
          Lsolver         = 'mumps'  #'ma97'  # Linear solver
          c_code          = False    # c_code
          shrinking_horizon = True
  
          return dt, x0, Lsolver, c_code, shrinking_horizon
  
      def DAE_system(self, uncertain_parameters=False):
          # Define vectors with names of states
          states     = ['x','n','q']
          nd         = len(states)
          xd         = SX.sym('xd',nd)
          for i in range(nd):
              globals()[states[i]] = xd[i]
  
          # Define vectors with names of algebraic variables
          algebraics = []
          na         = len(algebraics)
          xa         = SX.sym('xa',na)
          for i in range(na):
              globals()[algebraics[i]] = xa[i]
  
          set_points = []
          n_ref      = len(set_points)
          x_ref      = SX.sym('x_ref',n_ref)
          for i in range(n_ref):
              globals()[set_points[i]] = x_ref[i]
  
  
  
          # Define vectors with banes of input variables
          inputs     = ['L','Fn']
          nu         = len(inputs)
          u          = SX.sym("u",nu)
          for i in range(nu):
              globals()[inputs[i]] = u[i]
  
          # Define model parameter names and values
  
          #modparval = [0.0923*0.62, 178.85, 447.12, 393.10, 0.001, 504.49,
          #2.544*0.62*1e-4, 23.51, 800.0, 0.281, 16.89]
  
          if uncertain_parameters:
              modpar = ['u_m', 'k_s']
              nmp = len(modpar)
              theta = SX.sym("theta", nmp)
              for i in range(nmp):
                  globals()[modpar[i]] = theta[i]
          else:
              modpar = ['u_m', 'k_s']
              nmp = len(modpar)
              theta = SX.sym("theta", nmp)
              for i in range(nmp):
                  globals()[modpar[i]] = theta[i]
  
          k_i, K_N, u_d, Y_nx, k_m, k_sq, \
          k_iq, k_d, K_Np = 447.12, 393.10, 0.001, 504.49, \
                            2.544 * 0.62 * 1e-4, 23.51, 800.0, 0.281, 16.89
           # Additive measurement noise
      #    Sigma_v  = [400.,1e5,1e-2]*diag(np.ones(nd))*1e-6
  
          # Additive disturbance noise
      #    Sigma_w  = [400.,1e5,1e-2]*diag(np.ones(nd))*1e-6
  
          # Initial additive disturbance noise
      #    Sigma_w0 = [1.,150.**2,0.]*diag(np.ones(nd))*1e-3
          # Declare ODE equations (use notation as defined above)
  
          dx   = u_m * L/(L+k_s) * x * n/(n+K_N) - u_d*x
          dn   = - Y_nx*u_m* L/(L+k_s) * x * n/(n+K_N)+ Fn
          dq   = (k_m * L/(L+k_sq) * x - k_d * q/(n+K_Np))# * (sign(499.9999 - n)+1)/2 * (sign(x - 10.0001)+1)/2
  
                # * (tanh(1000*(500. - n))+1)/2 * (tanh(1000*(x - 10.000))+1)/2# * (sign(499.9999 - n)+1)/2 * (sign(x - 10.0001)+1)/2
  
          ODEeq =  [dx, dn, dq]
  
          # Declare algebraic equations
          Aeq = []
  
          # Define objective to be minimized
          t     = SX.sym('t')
          Obj_M = 0.#Function('mayer', [xd, x_ref, u], [-q])  # Mayer term
          Obj_L = 0.#Function('lagrange', [xd, x_ref, u], [0.])  # Lagrange term
          Obj_D = 0.#Function('Discrete', [xd, x_ref, u], [0.])  # Lagrange term
  
          R           = np.diagflat([3.125e-8, 3.125e-006])                         # Weighting of control penality
  
          # Define control bounds
          u_min = np.array([120., 0.])
          u_max = np.array([400., 40.])
          x_min = np.array([0.]*nd)
          x_max = np.array([np.inf]*nd)
          # Define constraint functions g(x) <= 0
          gequation = vertcat(n - 800.)#SX(0),SX(0))#
          ng = SX.size(gequation)[0]
          gfcn = Function('gfcn', [xd, xa, u], [gequation])
  
          return xd, xa, u, theta, ODEeq, Aeq, u_min, u_max, x_min, x_max, states,\
                 algebraics, inputs, nd, na, nu, n_ref, nmp, self.real_parameters, ng, gfcn, Obj_M, Obj_L, Obj_D, R
  
  
      def integrator_model(self):
          """
          This function constructs the integrator to be suitable with casadi environment, for the equations of the model
          and the objective function with variable time step.
           inputs: NaN
           outputs: F: Function([x, u, dt]--> [xf, obj])
          """
          dt, x0, Lsolver, c_code, shrinking_horizon = self.specifications()
  
          xd, xa, u, uncertainty, ODEeq, Aeq, u_min, u_max, x_min, x_max, states, \
          algebraics, inputs, nd, na, nu, n_ref, nmp, modparval, ng, gfcn, Obj_M, Obj_L, Obj_D, R = self.DAE_system()
  
          dae = {'x': vertcat(xd), 'z': vertcat(xa), 'p': vertcat(u, uncertainty),
                 'ode': vertcat(*ODEeq), 'alg': vertcat(*Aeq)}
          opts = {'tf': dt}  # interval length
          F = integrator('F', 'idas', dae, opts)
  
          return F
  
  