python - 尝试优化 Lugre 动态摩擦模型中的参数

标签 python numpy scipy physics ode

我收集了摩擦模型每个输出的 CSV 数据。该模型将表面之间的接触想象为一维刷毛,这些刷毛像 Spring 一样对弯曲使用react。摩擦力的模型为:

FL(V,Z) = sig0*Z +sig1*DZ/Dt +sig2*V 

其中 V 是表面速度 Z 是刷毛的偏转 DZ/Dt 是偏转率,等于:

DZ/Dt = V + abs(V)*Z/(Fc + (Fs-Fc)*exp(-(V^2/Vs^2))
      = V + abs(V)*Z/G(V)
      = V + H(V)*Z

其中 Fc 是运动中物体的摩擦力(常数),Fs 等于使物体运动所需的力(常数 > Fc),Vs 是在域之间转换所需的总速度(a我通过实验得出的常数)。 CSV 中提供了 block 的速度和位置以及摩擦力,所有这些都与时间有关。我还创建了一个易于积分的速度近似值作为时间的函数(三角函数)。

关于问题:代码与我尝试将列表传递给函数的方式不符(我认为)。

该函数传递的参数似乎可以工作(取自另一个仅绘制数据的文件),但是我尝试对 DZ/Dt 进行数值积分,并将 sig 参数拟合到导入的摩擦数据。

我导入的内容

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from scipy import optimize
import pylab as pp
from math import sin, pi, exp, fabs, pow

参数

Fc=2.7  #N
Fs=8.2  #N
Vs=.34  #mm/s

初始条件

ITime=Time[0]
Iz=[0,0,0]

建立摩擦模型

def velocity(time):
    V=-13/2*1/60*pi*sin(1/60*pi*time+pi)       
    return V

def g(v,vs,fc,fs,sig0):
    G=(1/sig0)*(fc+(fs-fc)*exp(-pow(v,2)/pow(vs,2)))
    return G

def h(v,vg):
  H=fabs(v)/vg
  return H

def findz(z, time, sig):
  Vx=velocity(time)
  VG=g(Vx,Vs,Fc,Fs,sig)
  HVx=h(Vx,VG)
  dzdt=Vx+HVx*z
  return dzdt

def friction(time,sig,iz):
  dz=lambda z,time: findz(z,time,sig)
  z=odeint(dz,iz,time)
  return sig[0]*z+sig[1]*findz(z,time,sig[0])+sig[2]*velocity(Time)

应该返回构造函数和数据之间的差异 生成包含优化参数的列表

def residual(sig):
  return Ff-friction(Time,sig,Iz)

SigG=[4,20,1]
SigVal=optimize.leastsq(residual,SigG)

print "parameter values are ",SigVal  

这会返回

line 56, in velocity
    V=-13/2*1/60*pi*sin(1/60*pi*time+pi)

TypeError: can't multiply sequence by non-int of type 'float'

这与我传递列表有关吗?

最佳答案

正如我在评论中提到的,Velocity() 是错误的原因,很可能是由于它使用时间值,而您传递了整个列表/数组(当您在 friction() 中调用它时,将其转换为 Velocity()(具有多个值)。

使用一些选定的值,并在缩短代码并传递 ITime 而不是 Time 后,代码可以正确运行,但由您来判断这是否是您分析的结果想要实现。下面是我的代码:

import numpy as np
from scipy import optimize
from scipy.integrate import odeint
from math import sin,  pi,  exp,  fabs

# Parameters
Fc = 2.7  #N
Fs = 8.2  #N
Vs = .34  #mm/s

# define test values for Ff and Time
Ff    = np.array([100, 50, 50])
Time  = np.array([10, 20, 30])

# Initial_conditions
ITime = Time[0]
Iz    = np.array([0, 0, 0])

# Building the friction model
V    = lambda                   t: (-13 / 2) * ( 1 / (60 * pi * sin(1 / 60 * pi * t + pi)))      
G    = lambda v, vs, fc, fs, sig0: (1 / sig0) * (fc + (fs - fc) * exp(-v**2 / vs**2))
H    = lambda               v, vg: fabs(v) / vg
dzdt = lambda         z,  t,  sig: V(t) + H(V(t), G(V(t), Vs, Fc, Fs, sig)) * z


def friction(t, sig, iz):
  dz = lambda z, t: dzdt(z, t, sig)
  z  = odeint(dz, iz, t)
  return sig[0]*z + sig[1]*dzdt(z, t, sig[0]) + sig[2]*V(t)

# Should return the difference between the Constructed function and the data
# and yield a list containing the optimized parameters

def residual(sig):

  return Ff-friction(ITime, sig, Iz)[0]

SigG   = np.array([4, 20, 1])
SigVal = optimize.leastsq(residual, SigG, full_output = False)

print("parameter values are ", SigVal )

输出:

parameter values are  (array([    4.        ,  3251.47271228, -2284.82881887]), 1)

关于python - 尝试优化 Lugre 动态摩擦模型中的参数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/56012137/

相关文章:

python - 在 google colab 上导入自己的 ipynb 文件

python - 将密集矩阵从文件直接读取到稀疏 numpy 数组中?

python - 如何使用 Python 插入 3D 曲面图(2D 数组)的缺失值(未定义区域)?

python - 在另一个 1d/2d 向量之后重采样(缩小)2D 向量

Python 数据矩阵检测

python - 何时以及为何使用 Django 开发服务器?

python - 使用(Python)NLTK提取产品核心关键词

python - 如何将 NumPy 数组元素从字符串更改为 int 或 float?

python - 由于numpy,Python脚本无法与Windows任务计划程序一起运行

python - 如何从数组中提取感兴趣的数字,同时将其与不同的数组进行比较?