model a logical statement function in Gekko

Question

I must create a friction model in GEKKO using the provided algorithm:

# u is the displacement variable
# vel is velocity , vel == u.dt()
# xr and fm are constant parameters


# update xs, frs:

if abs(vel) <= 1e-5:
    xs = u
    frs = f
else:
    xs = xs_prev
    frs = frs_prev

# update f

dx = u - xs
delta = frs/fm

if u == xs:
    f = frs
elif u > xs:
    f  =  frs + dx/(xr * (1-delta)+dx) * (fm - frs);
else:
    f  =  frs + dx/(xr * (1+delta)-dx) * (fm + frs);

I attempted various methods to simulate it in GEKKO but was unsuccess.

Approach 1: using Gekko logical functions

I added just the part of the code that deals with updating xs.

m = GEKKO()
m.time = np.linspace(0,10,100)
u = m.Var(0, fixed_initial=True)
t = m.Var(0, fixed_initial=True)
v = m.Var(3e-3 * 2*np.pi*0.1)

# previous time step 
xs_ = m.Var(0)
m.delay(xs, xs_, 1)

# Equations
m.Equation(t.dt()==1)
m.Equation(u == 3e-3 * m.sin(2*np.pi*0.1 * t))
m.Equation(v == u.dt())


# update xs
xs = m.Var(0, fixed_initial=True)
m.Equation(xs == m.if2(m.abs(v)<=1e-4, u,xs_) )

m.options.NODES = 2
m.options.IMODE = 6

m.solve(disp=False)

Using m.if2 and m.if3 raise following errors respectively:

 *** Error in syntax of function string: Missing operator
 
Position: 12                  
 v8-(abs(v3)<0.0001)
            ?

Exception:  @error: Model Expression
 *** Error in syntax of function string: Mismatched parenthesis
 
Position: 9                   
 (0.0001)))-((((1-int_v6))*(abs(v3))-slk_4
         ?

Approach 2: convert logical statements to smooth transition functions such as tanh and sigmoid.

This model sometimes works depending on the parameters.

m = GEKKO()
m.clear()
m.time = np.linspace(0,10,100)
u = m.Var(0)
up = m.Var(0)
t=m.Var(0)
m.Equation(t.dt()==1)
m.Equation(u == 3e-3 * m.sin(2*np.pi*0.1 * t))
xr = m.Param(0.2e-3)
xs = m.Var(0)
xs_ = m.Var(0)
fm = m.Param(value=250)
frs = m.Var(0)
f = m.Var(0)
frs_ = m.Var(0)
f_ = m.Var(0)
up_ = m.Var(0)
ds =m.Intermediate(u - xs)
dlt = m.Intermediate(frs/fm)

m.delay(xs, xs_, 1)
m.delay(frs, frs_, 1)
m.delay(f, f_, 1)
m.delay(up,up_,1)


#update xs and frs:
g_up = m.Intermediate(m.exp(-(up_/2e-4)**2))
m.Equation(u.dt()==up)

m.Equation(xs == u * g_up + xs_ * (1-g_up) )
m.Equation(frs == f_ * g_up + frs_ * (1-g_up) )
# 
# update f:
f_p = m.Intermediate(frs + ds/(xr * (1-dlt)+ds) * (fm - frs))
f_n = m.Intermediate(frs + ds/(xr * (1+dlt)-ds) * (fm + frs)) 
l0 = m.Intermediate(m.exp(-(ds/1e-9)**2))
l1 = m.Intermediate(m.sigmoid(1e6*ds))

m.Equation(f == l0*frs + l1*f_p + (1-l1)*f_n)

m.options.NODES = 2
m.options.IMODE = 6

m.solve(disp=False)

fig, ax = plt.subplots(3,1)
ax[0].plot(m.time, u, label='u')
ax[0].plot(m.time, xs, label='xs')
ax[0].legend()
ax[1].plot(m.time, f, label='f')
ax[1].plot(m.time, frs, label='f_')
ax[1].legend()
ax[2].plot(u, f, label='f')
ax[2].legend()

plt.show()

It doesn't work in all cases.

For example, if fm = 250, Gekko doesn't find a solution, but it works if fm = 300.

This is strange because the model does not consist of differential or Algebraic equations. It has to function with every parameters.

Here is the link for this model in Matlab Simulink.

What are the problems with both approaches?

Answer 1

Logical conditions in optimization with gradient-based optimizers need continuous 1st and 2nd derivatives. The sigmoid and exp functions that you developed also work, but the strong nonlinearity of these conditions:

l0 = m.Intermediate(m.exp(-(ds/1e-9)**2))
l1 = m.Intermediate(m.sigmoid(1e6*ds))

m.Equation(f == l0*frs + l1*f_p + (1-l1)*f_n)

means that the solver may fail to find a solution. Another way to model this is with Gekko logical conditions such as m.if2(), m.if3(), m.abs2(), and m.abs3(). Here is an example of a ramp on the force until the mass starts to move.

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

# Initialize the model
m = GEKKO(remote=True)

# Time horizon
n = 31
m.time = np.linspace(0,1,n)

# Parameters
mass = 5  # mass in kg
f_static = 1.0   # static friction coefficient
f_dynamic = 0.5  # dynamic friction coefficient

# Variables
v = m.Var(value=0)  # velocity, starting from rest

f_ramp = np.linspace(0,10,n)
force = m.Param(f_ramp)  # applied force

# Friction model: static and dynamic friction
epsilon = 1e-2  # small value to smooth the transition
fr = m.if2(m.abs2(v)-epsilon,f_static,f_dynamic)

# Add slack variable
u = m.Var(lb=0); slk = m.Var(lb=0)
m.Equation(u == force - fr * mass + slk)
m.Minimize(1e-3 * slk)

# Newton's second law with friction
m.Equation(mass * v.dt() == u)

# Solve the optimization problem
m.options.IMODE = 6 
m.options.SOLVER = 1
m.options.OTOL = 1e-8
m.options.RTOL = 1e-8
m.solve(disp=True)

# Plot the results
plt.figure(figsize=(6,4))
plt.subplot(2,1,1)
plt.plot(m.time, v.value, 'r--', label='Velocity')
plt.plot(m.time[1:], fr.value[1:], 'bo', label='Friction')
plt.ylabel('Velocity (m/s)'); plt.legend(); plt.grid()
plt.subplot(2,1,2)
plt.plot(m.time, force.value, 'k:', label='Force')
plt.ylabel('Force'); plt.legend(); plt.grid()
plt.xlabel('Time (s)')
plt.tight_layout()
plt.show()

Answer 2

I am new with gekko, but I think the usage of m.if2 might be wrong.

m.if2(m.abs(v)<=1e-4, u,xs_)

This should be written as:

m.if2(m.abs(v)-1e-4, u,xs_)

It means if m.abs(v)-1e-4 is less than 0, the value is u, or else it is xs_.