Defect in detect_simple_bounds?

[jgillis]

See Aachen example ex2_direct_collocation.py

[jgillis]

Does not occur for the pure SX counterpart

[jgillis]

The strategy to detect f2(p)x+f1(p)==0 structure from the 'simple' subset of g uses x=0 to evaluate. Non-simple entries of g get evaluated but their results (possibly NaN due to division by zero) are ignored. In other words, divisions by zero are localized. In FMU such errors are global.

Similar:

from casadi import *


opti = Opti()
x = opti.variable()

opti.minimize(x**2)

f = Function('f',[x],[3*x,x.attachAssert(x>1)**2],{"never_inline":True})

y,z = f(x)

opti.subject_to(y<=5)
opti.subject_to(z==3)

opti.set_initial(x, 3)

opti.solver("ipopt",{"detect_simple_bounds":True})

sol = opti.solve()

[jgillis]

In this particular case, simplification of vertcat(stuff,stuff,stuff,daefun())[:3] to vertcat(stuff,stuff,stuff) would work.

## 2.1
import casadi as ca
import numpy as np
import os
from fmi_helpers import unpack_fmu, cleanup_fmu

# Unpack the FMU
unpacked_fmu = unpack_fmu('resources/bvsm.fmu')

# Load it into CasADi, display
dae = ca.DaeBuilder('BVSM', unpacked_fmu)
dae.disp(True)

# Make k1 and k2 parameters instead of controls
dae.set_all('u', ['Qf'])
dae.set_all('p', ['k1', 'k2'])
dae.disp(True)

## 2.2

# Create a function for evaluating the system dynamics
daefun = dae.create('daefun')
print(daefun)

# Also create an output function
yfun = dae.create('yfun', ['x','u', 'p'], ['ydef'])
print(yfun)

# Evaluate the output function at some point
y_test = yfun([2000., 2., 0.02, 0.2, 1.5], 0.5, [1000., 1000.])
print(y_test)

## 2.3

# Load actuator and measurement data
from scipy.io import loadmat
data = loadmat('resources/bvsm.mat')
tout = data['tout'].flatten()
lc_m = data['lc_m']
Qf = data['Qf']
N = len(tout) - 1

# Create Integrator instance, simulate the system
simulator = ca.integrator('simulator', 'cvodes', daefun, 0, tout)
x0 = dae.start(dae.x())  # Initial state
p0 = dae.start(dae.p())  # Parameter guess
u0 = Qf                  # Control (given)
res = simulator(x0 = x0, p = p0, u = u0)
x_sim = res['xf']

# Evaluate the outputs, plot
y_all = yfun(x_sim, u0, p0)


## 2.4

# Recreate a new simulator over 10 min
F = ca.integrator('F', 'collocation', daefun, 0, 10,
    dict(number_of_finite_elements=1))

# Simulate the system using this integrator instead

x_sim = [x0]
xk = x0
for k in range(90):
  Fk = F(x0 = xk, p = p0, u = u0[:, k])
  xk = Fk['xf']
  x_sim.append(xk)
x_sim = ca.hcat(x_sim)

# Evaluate the outputs
y_all = yfun(x_sim, u0, p0)

## 2.5
# Start with an empty NLP
opti = ca.Opti()



# Create a decision variable for nA0
# Allow the value to vary between 80 % and 120 % of the given value
nA0 = opti.variable()
nA_guess = dae.start('nA')
nA_guess = 2.416485390731904
opti.set_initial(nA0, nA_guess)
opti.subject_to(0.8*nA_guess <= (nA0 <= 1.2*nA_guess))
# Symbolic initial conditions
xk = ca.MX(x0)
xk[dae.x().index('nA')] = nA0

# Formulate the NLP with single shooting
J = 0
y_sim = []
x_sim = [xk]
for k in range(N):
    # Get control input
    uk = u0[:, k]
    # Call the integrator symbolically
    Fk = F(x0 = xk, p = p0, u = uk)
    xk = Fk['xf']
    # Evaluate output function
    yk = yfun(xk, uk, p0)
    # Add objective terms
    lc_k = lc_m[:, k + 1]
    if not np.isnan(lc_k): J = J + ((lc_k - yk) / yk) ** 2
    y_sim.append(yk)
    x_sim.append(xk)
y_sim = ca.hcat(y_sim)
x_sim = ca.hcat(x_sim)

# Optimize
opti.minimize(J)
opti.solver('ipopt')
sol = opti.solve()

x_sim = sol.value(x_sim)
print(sol.value(x_sim[:,-1]))

# Get the optimal parameters
nA0_opt = sol.value(nA0)
print(f'Solution: nA0 = {nA0_opt}')

# Evaluate the outputs at the optimal solution
y_all = sol.value(y_sim)

## 2.5 with k1 and k2 free

# Create a new NLP with k1 and k2 free
opti = ca.Opti()

# Create a decision variable for nA0
# Allow the value to vary between 80 % and 120 % of the given value
nA0 = opti.variable()
nA_guess = dae.start('nA')
opti.set_initial(nA0, nA_guess)
opti.subject_to(0.8*nA_guess <= (nA0 <= 1.2*nA_guess))

# Create decision variable for k1
k1 = opti.variable()                 # New decision variable
opti.set_initial(k1, dae.start('k1')) # Initial guess
opti.subject_to(dae.min('k1') <= (k1 <= dae.max('k1')))

# Create decision variable for k2
k2 = opti.variable()                 # New decision variable
opti.set_initial(k2, dae.start('k2')) # Initial guess
opti.subject_to(dae.min('k2') <= (k2 <= dae.max('k2')))

# Symbolic parameter vector
p = ca.vertcat(k1, k2)

degree = 3
method = 'radau'
dt = 10
tau = ca.collocation_points(degree,method)
[C,D,B] = ca.collocation_coeff(tau)

# Symbolic initial conditions
xk = ca.MX(x0)
xk[dae.x().index('nA')] = nA0

# Formulate the NLP with single shooting
J = 0
y_sim = []
for k in range(N):
    # Get control input
    uk = u0[:, k]

    # Decision variables for helper states at each collocation point
    Xc = opti.variable(5, degree)
    Z = ca.horzcat(xk,Xc)
    Pidot = (Z @ C)/dt
    opti.subject_to(Pidot==daefun(x=Xc,u=uk,p=p)["ode"])
    
    # Continuity constraints
    xk_next = opti.variable(5)
    opti.subject_to(Z @ D==xk_next)
    opti.set_initial(Xc, ca.repmat(x_sim[:,k],1,degree))
    opti.set_initial(xk_next, x_sim[:,k+1])

    # Evaluate output function
    yk = yfun(xk_next, uk, p)
    # Add objective terms
    lc_k = lc_m[:, k + 1]
    if not np.isnan(lc_k):
        J = J + ((lc_k - yk) / yk) ** 2
        
    xk = xk_next
    y_sim.append(yk)
y_sim = ca.hcat(y_sim)

# Optimize
opti.minimize(J)
opti.solver('ipopt',{"detect_simple_bounds":True,"common_options.final_options":{"print_in":True}})
sol = opti.solve()

# Get the optimal parameters
k1_opt = sol.value(k1)
k2_opt = sol.value(k2)
nA0_opt = sol.value(nA0)
print(f'Solution: k1 = {k1_opt}, k2 = {k2_opt}, nA0 = {nA0_opt}')

# Evaluate the outputs at the optimal solution
y_all = sol.value(y_sim)