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)
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.

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"