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

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_.