我正在 Python 中使用 GEKKO 来估计弹跳球的轨迹。为此，我需要估计 2 个变量：e_1（恢复系数）和 q_1（每次弹跳时的水平速度损失）。
我已为其编写了以下代码，但参数似乎没有更新，尽管解算器已成功执行。参数的初始值与参数的最终优化值相同
e_1= 0.8
和
q_1= 1
代码：
import numpy as np
from gekko import GEKKO
# Actual trajectories
# actually y
act_x = np.array([0.019053200000000013, 0.08770320000000008, 0.1850550000000004, 0.2825550000000008, 0.38005500000000114, 0.4775550000000015, 0.5750549999999806, 0.6725549999999532, 0.7700549999999258, 0.8675549999998984, 0.965054999999871, 1.0625549999998436, 1.1600549999998162, 1.2575549999997888, 1.3550549999997614, 1.452554999999734, 1.5500549999997066, 1.6475549999996792, 1.7450549999996519, 1.8425549999996245, 1.940054999999597, 2.0375549999996125, 2.1234107142854164, 2.193053571428365, 2.2626964285713136, 2.332339285714262, 2.401982142857211, 2.4716250000001594, 2.541267857143108, 2.6109107142860566, 2.680553571429005, 2.750196428571954, 2.8198392857149024, 2.889482142857851, 2.9591250000007996, 3.028767857143748, 3.098410714286697, 3.1680535714296454, 3.237696428572594, 3.3073392857155426, 3.376982142858491, 3.44662500000144, 3.5162678571443884, 3.585910714287337, 3.6555535714302856, 3.725196428573234, 3.794839285716183, 3.8644821428591314, 3.93412500000208, 4.003767857145029, 4.073410714287977, 4.143053571430926, 4.212696428573874, 4.282339285716823, 4.351982142859772, 4.42162500000272, 4.491267857145669, 4.560910714288617, 4.630553571431566, 4.700196428574515, 4.769839285717463, 4.839482142860412, 4.90912500000336, 4.978767857146309, 5.048410714289258, 5.118053571432206, 5.187696428575155, 5.257339285718103, 5.326982142861052, 5.396625000004001, 5.466267857146949, 5.535910714289898, 5.605553571432846, 5.675196428575795, 5.744839285718744, 5.814482142861692, 5.884125000004641, 5.953767857147589, 6.023410714290538, 6.093053571433487, 6.162696428576435, 6.232339285719384, 6.301982142862332, 6.371625000005281, 6.44126785714823, 6.510910714291178, 6.580553571434127, 6.650196428577075, 6.719839285720024, 6.789482142862973, 6.859125000005921, 6.92876785714887, 6.998410714291818, 7.068053571434767, 7.137696428577716, 7.207339285720664, 7.276982142863613, 7.346625000006561, 7.41626785714951, 7.485910714292459, 7.555553571435407, 7.625196428578356, 7.694839285721304, 7.764482142864253])
# actually z
act_y = np.array([1.0379079699999998, 1.1754534700000003, 1.3602534700000006, 1.5209239699999955, 1.6570944699999859, 1.7687649699999706, 1.8559354699999533, 1.9186059699999365, 1.95677646999992, 1.970446969999904, 1.9596174699998887, 1.9242879699998736, 1.8644584699998585, 1.7801289699998444, 1.6712994699998283, 1.537969969999807, 1.3801404699997812, 1.19781096999975, 0.9909814699997134, 0.7596519699996719, 0.5038224699996257, 0.22349296999957416, 0.13930352847394978, 0.3499050645343172, 0.5360066005946793, 0.6976081366550367, 0.8347096727153887, 0.9473112087757359, 1.0354127448360804, 1.0990142808964258, 1.1381158169567716, 1.1527173530171173, 1.1428188890774638, 1.1084204251378103, 1.0495219611981568, 0.9661234972585043, 0.8582250333188504, 0.7258265693791917, 0.5689281054395274, 0.3875296414998584, 0.18163117756018424, 0.11887053865291446, 0.2798887215524012, 0.41640690445188283, 0.5284250873513592, 0.6159432702508338, 0.6789614531503088, 0.7174796360497837, 0.7314978189492595, 0.7210160018487353, 0.6860341847482112, 0.6265523676476875, 0.5425705505471646, 0.4340887334466397, 0.3011069163461103, 0.14362509924557573, 0.10758844178113687, 0.2208353837910192, 0.309582325800899, 0.37382926781077935, 0.41357620982065996, 0.42882315183054115, 0.4195700938404224, 0.385817035850304, 0.3275639778601862, 0.24481091987006875, 0.13755786187995, 0.07905142871570557, 0.1646809455657065, 0.22581046241570793, 0.2624399792657099, 0.2745694961157122, 0.26219901296571463, 0.2253285298157173, 0.16395804666572056, 0.07808756351572432, 0.09996825173824968, 0.1511630228129375, 0.17785779388762588, 0.1800525649623144, 0.15774733603700322, 0.11094210711169239, 0.05563702598852237, 0.09995862508974714, 0.11978022419097237, 0.1151018232921977, 0.08592342239342333, 0.06022248147389577, 0.08484848858531084, 0.0849744956967261, 0.06060050280814157, 0.06605277311641561, 0.06792022730639775, 0.05243051383233173, 0.060240252169206004, 0.0532114047328211, 0.05285089806699147, 0.052793140813351076, 0.05130362411928389, 0.05057861286791727, 0.049998536419320685, 0.049997973905616416, 0.049997920093917625, 0.049997920094010155])
def get_initial_trajectory(e_1, q_1):
# Constants
g = 9.81 # gravity (m/s^2)
x0, y0 = 0, 1 # initial position (m)
vx0 = 1.95 # Velocity in x direction
vy0 = 4.4 # Velocity in y direction
t_step = 0.01 # time step for simulation (s)
t_max = 5 # max time for simulation (s)
x_vals = [x0]
y_vals = [y0]
t = 0
x, y = x0, y0
vx, vy = vx0, vy0
while t < t_max:
t += t_step
x += vx * t_step
y += vy * t_step - 0.5 * g * t_step ** 2
vy -= g * t_step
if y <= 0: # Ball hits the ground
y = 0
vy = -e_1 * vy
vx = vx * q_1
x_vals.append(x)
y_vals.append(y)
if len(x_vals) > 1 and x_vals[-1] == x_vals[-2] and y_vals[-1] == y_vals[-2]:
break
return np.array(x_vals, dtype='float64'), np.array(y_vals, dtype='float64')
def cost_function(e_1, q_1):
x_vals, y_vals = get_initial_trajectory(e_1, q_1)
y_temp = np.interp(act_x, x_vals, y_vals)
return np.sum((act_y - y_temp)**2)
# Optimization with Gekko
m = GEKKO(remote=False)
# Define variables
e_1 = m.Var(value=0.8, name='e_1')
q_1 = m.Var(value=1, name='q_1')
m.Minimize(cost_function(e_1.value, q_1.value))
# Set the objective
m.options.SOLVER=3
# Solve the optimization problem
m.solve(disp=True)
# Print results
print(f"Estimated e_1: {e_1.value[0]}")
print(f"Estimated q_1: {q_1.value[0]}")
输出：
----------------------------------------------------------------
APMonitor, Version 1.0.3
APMonitor Optimization Suite
----------------------------------------------------------------
--------- APM Model Size ------------
Each time step contains
Objects : 0
Constants : 0
Variables : 2
Intermediates: 0
Connections : 0
Equations : 1
Residuals : 1
________________________________________________
WARNING: objective equation 1 has no variables
ss.Eqn(1)
0 = 12.347769475360362
________________________________________________
Number of state variables: 2
Number of total equations: - 0
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 2
solver 3 not supported
using default solver: APOPT
----------------------------------------------
Steady State Optimization with APOPT Solver
----------------------------------------------
1 1.23478E+01 0.00000E+00
Successful solution
---------------------------------------------------
Solver : IPOPT (v3.12)
Solution time : 2.689999999711290E-002 sec
Objective : 12.3477694753604
Successful solution
---------------------------------------------------
Estimated e_1: 0.8
Estimated q_1: 1.0
Process finished with exit code 0
如果需要任何其他信息，请告诉我。
我发现了代码中导致参数没有更新的几个问题：
-
传递给
m.Minimize()的目标函数不正确。 你应该传递一个 Gekko 表达式作为目标函数，而不是传递一个已经计算好的数值。当前，cost_function(e_1.value, q_1.value)在每次迭代中只计算一次，并返回一个常数，因此 Gekko 无法对其进行优化。 -
get_initial_trajectory()函数中使用了 Gekko 变量，但没有将其作为 Gekko 表达式处理。 你需要使用 Gekko 的运算符来执行计算，以便 Gekko 可以跟踪变量之间的关系并计算导数。
以下是修改后的代码：
import numpy as np
from gekko import GEKKO
# ... (act_x 和 act_y 数据保持不变) ...
def get_initial_trajectory(m, e_1, q_1):
# ... (常数定义保持不变) ...
x_vals = [x0]
y_vals = [y0]
t = 0
x = m.Var(value=x0)
y = m.Var(value=y0)
vx = m.Var(value=vx0)
vy = m.Var(value=vy0)
m.Equation(x.dt() == vx)
m.Equation(y.dt() == vy - 0.5 * g * t_step)
m.Equation(vy.dt() == -g)
for _ in range(int(t_max / t_step)):
m.time.push(t)
t += t_step
m.solve(disp=False)
# 使用 Gekko 的 if3 函数来模拟条件语句
vy = m.if3(y - 1e-6, vy, -e_1 * vy)
vx = m.if3(y - 1e-6, vx, vx * q_1)
# 更新 x 和 y 的值
x_vals.append(x.value[0])
y_vals.append(y.value[0])
return np.array(x_vals), np.array(y_vals)
def cost_function(m, e_1, q_1):
x_vals, y_vals = get_initial_trajectory(m, e_1, q_1)
y_temp = np.interp(act_x, x_vals, y_vals)
return np.sum((act_y - y_temp)**2)
# Optimization with Gekko
m = GEKKO(remote=False)
# Define variables
e_1 = m.Var(value=0.8, name='e_1', lb=0.1, ub=1)
q_1 = m.Var(value=1, name='q_1', lb=0.1, ub=1)
# 计算初始轨迹
x_init, y_init = get_initial_trajectory(m, e_1, q_1)
# 使用 Gekko 的数组函数来定义目标函数
m.Minimize(sum([(act_y[i] - y_init[i])**2 for i in range(len(act_x))]))
# Set the objective
m.options.IMODE = 6
m.options.SOLVER=1
# Solve the optimization problem
m.solve(disp=True)
# Print results
print(f"Estimated e_1: {e_1.value[0]}")
print(f"Estimated q_1: {q_1.value[0]}")
主要改动:
-
将 Gekko 模型
m作为参数传递给get_initial_trajectory()和cost_function()。 -
在
get_initial_trajectory()中，使用 Gekko 变量和方程式来计算轨迹，而不是使用 Python 的基本运算符。 -
使用
m.if3()函数来模拟弹跳时的条件语句，这对于 Gekko 的求解器是必要的。 - 使用 Gekko 的数组函数来计算目标函数。
-
将
m.options.IMODE设置为6以进行参数估计。
这些修改应该可以解决你的问题，并使 Gekko 能够正确地更新参数
e_1
和
q_1
。