text: # GP NMPC problem setup
  import numpy as np
  from casadi import *
  
  
  class model1:
      def __init__(self):
          self.real_parameters = [1., 1., 1.]
  
      def specifications(self):
          ''' Specify Problem parameters '''
          tf = 2.5  # final time
          nk = 10  # sampling points
          dt = tf / nk
          x0 = np.array([0.95, 0.04])
          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     = ['x1','x2']
          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     = ['u']
          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 = ['p1', 'p2', 'p3']
              nmp = len(modpar)
              theta = SX.sym("theta", nmp)
              for i in range(nmp):
                  globals()[modpar[i]] = theta[i]
          else:
              modpar = ['p1','p2', 'p3']
              nmp = len(modpar)
              theta = SX.sym("theta", nmp)
              for i in range(nmp):
                  globals()[modpar[i]] = theta[i]
  
          dx1  = (p1 + p3)*(x2-1)*x1 + (p2+u)*x2
          dx2  = p1*(1-x2)*x1-(p2+u)*x2
  
  
          ODEeq =  [dx1, dx2]
  
          # 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([0.])
          u_max = np.array([5.])
          x_min = np.array([0.]*nd)
          x_max = np.array([np.inf]*nd)
          # Define constraint functions g(x) <= 0
          gequation = vertcat(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
  
  
  class model2:
      def __init__(self):
          self.real_parameters =  [1., 1., 1.]
      def specifications(self):
          ''' Specify Problem parameters '''
          tf              = 2.5      # final time
          nk              = 10        # sampling points
          dt              = tf/nk
          x0              = np.array([0.95,0.04])
          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     = ['x1','x2']
          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     = ['u']
          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 = ['p1', 'p2', 'p3']
              nmp = len(modpar)
              theta = SX.sym("theta", nmp)
              for i in range(nmp):
                  globals()[modpar[i]] = theta[i]
          else:
              modpar = ['p1', 'p2', 'p3']
              nmp = len(modpar)
              theta = SX.sym("theta", nmp)
              for i in range(nmp):
                  globals()[modpar[i]] = theta[i]
  
  
          dx1   = p1*(x2-1)*x1 + (p2+u)*x2
          dx2   = p1*(1-x2)*x1-((p2+u)+p3)*x2
  
  
  
          ODEeq =  [dx1, dx2]
  
          # 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([0.])
          u_max = np.array([5.])
          x_min = np.array([0]*nd)
          x_max = np.array([np.inf]*nd)
          # Define constraint functions g(x) <= 0
          gequation = vertcat(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
  
  