Jump to content

Three Tank multimode problem (python/casadi)

From mintOC
Revision as of 09:17, 14 March 2020 by ClemensZeile (talk | contribs) (Created page with "=== Casadi === The model in Python code for a fixed control discretization grid using direct multiple shooting. We use pycombina for solving the (CIA) problem. <source lang="...")
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

Casadi

The model in Python code for a fixed control discretization grid using direct multiple shooting. We use pycombina for solving the (CIA) problem. <source lang="Python">

  1. Three tank optimal control
  2. ----------------------
  3. An optimal control problem (OCP),
  4. solved with direct multiple-shooting.
  6. For more information on CasADi see: http://labs.casadi.org/OCP

from casadi import * import pylab as pl from pycombina import BinApprox, CombinaBnB import matplotlib.pyplot as plt

                          1. Input variables #############

T = 12 N = 100 # number of control intervals


                          1. End of Input variables #############

opti = Opti() # Optimization problem


  1. ------- parameter values --------------
  2. Model parameters

k_1 = 2 k_2 = 3 k_3 = 1 k_4 = 3


c_1 = 1 c_2 = 2 c_3 = 0.8 T = 12

  1. ---- decision variables ---------

X = opti.variable(4,N+1) # state trajectory x1 = X[0,:] x2 = X[1,:] x3 = X[2,:] x4 = X[3,:]


U = opti.variable(3,N) # control trajectory (throttle) u1 = U[0,:] u2 = U[1,:] u3 = U[2,:]

  1. ---- objective ---------

opti.minimize(x4[N]) # Objective as Mayer term

  1. ---- dynamic constraints -------

f = lambda x,u: vertcat(c_1*u[0]+c_2*u[1]-(x[0])**0.5- u[2]* (c_3*x[0])**0.5 , (x[0])**0.5- (x[1])**0.5, (x[1])**0.5- (x[2])**0.5 +u[2]* (c_3*x[0])**0.5 , u[1]*0.0+ k_1*(x[1]-k_2)**2 + k_3*(x[2]-k_4)**2 ) # dx/dt = f(x,u)


dt = T/N # length of a control interval for k in range(N): # loop over control intervals 

  # Runge-Kutta 4 integration
  k1 = f(X[:,k],         U[:,k])
  k2 = f(X[:,k]+dt/2*k1, U[:,k])
  k3 = f(X[:,k]+dt/2*k2, U[:,k])
  k4 = f(X[:,k]+dt*k3,   U[:,k])
  x_next = X[:,k] + dt/6*(k1+2*k2+2*k3+k4) 
  opti.subject_to(X[:,k+1]==x_next) # close the gaps
  1. ---- path constraints -----------

opti.subject_to(opti.bounded(0,U,1)) # control is limited opti.subject_to(u1+u2+u3 == 1) # SOS1 constraint


  1. ---- boundary conditions --------

opti.subject_to(x1[0]==2) # intitial values opti.subject_to(x2[0]==2) # ... opti.subject_to(x3[0]==2) # ... opti.subject_to(x4[0]==0) # ... mayer term auxiliary state

  1. opti.subject_to(pos[-1]==1) # finish line at position 1
  1. ---- misc. constraints ----------
  2. opti.subject_to(T>=0) # Time must be positive
  1. ---- initial values for solver ---

opti.set_initial(x1, 2) opti.set_initial(x2, 2) opti.set_initial(x3, 2) opti.set_initial(x4, 0) opti.set_initial(u1, 1) opti.set_initial(u2, 0) opti.set_initial(u3, 0)


  1. ---- solve NLP ------

opti.solver("ipopt") # set numerical backend sol = opti.solve() # actual solve


  1. ---- post-processing ------

from pylab import plot, step, figure, legend, show, spy

plt.plot(pl.linspace(0, 12, num=N+1),sol.value(x1),label="x1") plt.plot(pl.linspace(0, 12, num=N+1),sol.value(x2),label="x2") plt.plot(pl.linspace(0, 12, num=N+1),sol.value(x3),label="x3")

  1. plot(limit(sol.value(pos)),'r--',label="speed limit")

plt.step(pl.linspace(0, 12, num=N),sol.value(U[0,:]),'k',label="a1") plt.step(pl.linspace(0, 12, num=N),sol.value(U[1,:]),label="a2") plt.step(pl.linspace(0, 12, num=N),sol.value(U[2,:]),label="a3")

plt.xlabel('t') plt.ylabel('state and control values') plt.title('Differential state trajectories for relaxed controls')

  1. xlabel=("t")

plt.legend(loc="upper left")

  1. title="

plt.savefig('three_tank_relaxed_solution.png')

      1. Setting up CIA problem
################################################

t = np.linspace(0,T,N+1) b_rel = np.array([[min(1,max(0,x)) for x in sol.value(U[0,:])], [min(1,max(0,x)) for x in sol.value(U[1,:])], [min(1,max(0,x)) for x in sol.value(U[2,:])]])


binapprox = BinApprox(t = t, b_rel = b_rel, binary_threshold = 1e-3, \ 

   off_state_included = False)
  1. binapprox.set_n_max_switches(n_max_switches = max_switches)
  2. binapprox.set_min_up_times(min_up_times = [min_up_time])
  3. binapprox.set_min_down_times(min_down_times = [min_down_time])

combina = CombinaBnB(binapprox) combina.solve(use_warm_start=False, bnb_search_strategy='dfs', bnb_opts={'max_iter': 1000000})

b_bin_orig = pl.asarray(binapprox.b_bin)


      1. Re-run Optimization Problem with fixed binary controls
      2. 'Dirty hack': just copy previous problem set up due to reuse issues of constraints

sol1_x1 = sol.value(x1) sol1_x2 = sol.value(x2) sol1_x3 = sol.value(x3)


opti2 = Opti() # Optimization problem

  1. ---- decision variables ---------

X = opti2.variable(4,N+1) # state trajectory x1 = X[0,:] x2 = X[1,:] x3 = X[2,:] x4 = X[3,:]

U = opti2.variable(3,N) # control trajectory (throttle) u1 = U[0,:] u2 = U[1,:] u3 = U[2,:]

  1. ---- objective ---------

opti2.minimize(x4[N]) # Objective as Mayer term

  1. ---- dynamic constraints --------

f = lambda x,u: vertcat(c_1*u[0]+c_2*u[1]-(x[0])**0.5- u[2]* (c_3*x[0])**0.5 , (x[0])**0.5- (x[1])**0.5, (x[1])**0.5- (x[2])**0.5 +u[2]* (c_3*x[0])**0.5 , u[1]*0.0+ k_1*(x[1]-k_2)**2 + k_3*(x[2]-k_4)**2)


dt = T/N # length of a control interval for k in range(N): # loop over control intervals 

  # Runge-Kutta 4 integration
  k1 = f(X[:,k],         U[:,k])
  k2 = f(X[:,k]+dt/2*k1, U[:,k])
  k3 = f(X[:,k]+dt/2*k2, U[:,k])
  k4 = f(X[:,k]+dt*k3,   U[:,k])
  x_next = X[:,k] + dt/6*(k1+2*k2+2*k3+k4) 
  opti2.subject_to(X[:,k+1]==x_next) # close the gaps
  1. ---- path constraints -----------

opti2.subject_to(U==b_bin_orig) # control is fixed

  1. ---- boundary conditions --------

opti2.subject_to(x1[0]==2) # intitial values opti2.subject_to(x2[0]==2) # ... opti2.subject_to(x3[0]==2) # ... opti2.subject_to(x4[0]==0) # ... mayer term auxiliary state

  1. opti.subject_to(pos[-1]==1) # finish line at position 1
  1. ---- misc. constraints ----------
  2. opti.subject_to(T>=0) # Time must be positive
  1. ---- initial values for solver ---

opti2.set_initial(sol.value_variables())

  1. ---- initial values for solver ---

opti2.set_initial(x1, sol1_x1) opti2.set_initial(x2, sol1_x2) opti2.set_initial(x3, sol1_x3) opti2.set_initial(x4, 0)

  1. opti.set_initial(u1, 0)
  2. opti.set_initial(u2, 0)
  3. opti.set_initial(u3, 1)


  1. ---- solve NLP ------

opti2.solver("ipopt") # set numerical backend sol2 = opti2.solve() # actual solve


        2. Plot the binary feasible trajectory solution

figure()

plt.plot(pl.linspace(0, 12, num=N+1),sol2.value(x1),label="x1") plt.plot(pl.linspace(0, 12, num=N+1),sol2.value(x2),label="x2") plt.plot(pl.linspace(0, 12, num=N+1),sol2.value(x3),label="x3")

  1. plot(limit(sol.value(pos)),'r--',label="speed limit")

plt.step(pl.linspace(0, 12, num=N),sol2.value(U[0,:]),'k',label="w1") plt.step(pl.linspace(0, 12, num=N),sol2.value(U[1,:]),label="w2") plt.step(pl.linspace(0, 12, num=N),sol2.value(U[2,:]),label="w3")

plt.xlabel('t') plt.ylabel('state and control values') plt.title('Differential state trajectories for binary controls')

  1. xlabel=("t")

plt.legend(loc="upper left")

  1. title="

show() plt.savefig('three_tank_binary_solution.png')

        2. Plot also the (CIA) rounding trajectories

f, (ax1, ax2, ax3) = pl.subplots(3, sharex = True)

ax1.step(t[:-1], b_bin_orig[0,:], label = "w", color = "C0", where = "post")

  1. ax1.step(t[:-1], b_bin_red[0,:], label = "b_bin_red", color = "C0", linestyle = "dashed", where = "post")

ax1.scatter(t[:-1], b_rel[0,:], label = "a", color = "C0", marker = "x") ax1.legend(loc = "upper left") ax1.set_ylabel("w_1")

ax2.step(t[:-1], b_bin_orig[1,:], label = "w", color = "C1", where = "post")

  1. ax2.step(t[:-1], b_bin_red[1,:], label = "b_bin_red", color = "C1", linestyle = "dashed", where = "post")

ax2.scatter(t[:-1], b_rel[1,:], label = "a", color = "C1", marker = "x") ax2.legend(loc = "upper left") ax2.set_ylabel("w_2")

ax3.step(t[:-1], b_bin_orig[2,:], label = "w", color = "C2", where = "post")

  1. ax3.step(t[:-1], b_bin_red[2,:], label = "b_bin_red", color = "C2", linestyle = "dashed", where = "post")

ax3.scatter(t[:-1], b_rel[2,:], label = "a", color = "C2", marker = "x") ax3.legend(loc = "upper left") ax3.set_ylabel("w_3") ax3.set_xlabel("t")

show() plt.savefig('three_tank_rounding_solution.png') 

################################################


<source>

Retrieved from "https://mintoc.de/index.php?title=Three_Tank_multimode_problem_(python/casadi)&oldid=2335"