[Python] model a logical statement function in Gekko

I must create a friction model 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?

Continue reading...