问题描述
高尔夫击球是一个经典的物理优化问题：球员需要选择击球的初始速度和发射角度，使球落入球洞中。球在空中受到重力和空气阻力的影响，轨迹呈抛物线（或受阻力修正的抛物线）。
本文将用Java实现高尔夫击球的物理模拟系统，并通过梯度下降算法自动寻找最优击球角度和力度，使球尽可能准确地落入目标区域。
核心算法思路
物理模型
在不考虑空气阻力的情况下，球的轨迹满足抛体运动方程：
$$x(t) = v_0 \cos(\theta) \cdot t$$
$$y(t) = v_0 \sin(\theta) \cdot t – \frac{1}{2} g t^2$$
其中 $v_0$ 为初速度，$\theta$ 为发射角度，$g = 9.8 \text{m/s}^2$ 为重力加速度。
考虑空气阻力时，阻力与速度平方成正比：
$$F_d = -k \cdot v^2 \cdot \frac{\vec{v}}{|v|}$$
梯度下降优化
目标函数为球落点与球洞位置的欧氏距离：
$$L(v_0, \theta) = \sqrt{(x_{land} – x_{hole})^2 + (y_{land} – y_{hole})^2}$$
通过计算目标函数对 $v_0$ 和 $\theta$ 的偏导数，迭代更新参数：
$$v_0 \leftarrow v_0 – \alpha \frac{\partial L}{\partial v_0}$$
$$\theta \leftarrow \theta – \alpha \frac{\partial L}{\partial \theta}$$
由于解析导数难以求得，本文采用数值差分近似梯度。
Java实现
import java.util.*;
/**
* 高尔夫击球模拟与最优策略求解
* 使用抛物线物理模型和梯度下降寻找最优击球参数
*/
public class GolfSimulator {
// 物理常量
private static final double G = 9.8; // 重力加速度 (m/s^2)
private static final double AIR_RESISTANCE = 0.001; // 空气阻力系数
private static final double TIME_STEP = 0.01; // 仿真时间步长 (s)
private static final double GROUND_Y = 0.0; // 地面高度
// 球场参数
private double holeX; // 球洞水平距离 (m)
private double holeY; // 球洞深度/高度 (m)，通常在地面上为0
private double holeRadius; // 球洞半径 (m)
public GolfSimulator(double holeX, double holeY, double holeRadius) {
this.holeX = holeX;
this.holeY = holeY;
this.holeRadius = holeRadius;
}
/**
* 球的状态向量 [x, y, vx, vy]
*/
public static class BallState {
public double x, y, vx, vy;
public BallState(double x, double y, double vx, double vy) {
this.x = x; this.y = y; this.vx = vx; this.vy = vy;
}
public BallState copy() {
return new BallState(x, y, vx, vy);
}
}
/**
* 无阻力抛物线运动的解析解
* @param v0 初速度
* @param angle 发射角度（弧度）
* @return 落点水平距离
*/
public double analyticalRange(double v0, double angle) {
// 飞行时间 t = 2 * v0 * sin(theta) / g
double flightTime = 2 * v0 * Math.sin(angle) / G;
// 射程 R = v0^2 * sin(2*theta) / g
return v0 * Math.cos(angle) * flightTime;
}
/**
* 带空气阻力的数值仿真
* @param v0 初速度
* @param angle 发射角度（弧度）
* @return 落点状态
*/
public BallState simulateWithDrag(double v0, double angle) {
BallState state = new BallState(0, 0.01, // 略高于地面避免数值问题
v0 * Math.cos(angle), v0 * Math.sin(angle));
while (state.y > GROUND_Y || state.vy > 0) {
// 计算速度大小
double v = Math.sqrt(state.vx * state.vx + state.vy * state.vy);
// 阻力加速度
double dragAx = -AIR_RESISTANCE * v * state.vx;
double dragAy = -AIR_RESISTANCE * v * state.vy;
// 更新速度
state.vx += dragAx * TIME_STEP;
state.vy += (dragAy - G) * TIME_STEP;
// 更新位置
state.x += state.vx * TIME_STEP;
state.y += state.vy * TIME_STEP;
}
// 线性插值精确落点
if (state.y < GROUND_Y) {
double ratio = (state.y - GROUND_Y) / state.vy;
state.x -= state.vx * ratio;
state.y = GROUND_Y;
}
return state;
}
/**
* 目标函数：落点与球洞的距离
*/
public double lossFunction(double v0, double angle) {
BallState landing = simulateWithDrag(v0, angle);
double dx = landing.x - holeX;
double dy = landing.y - holeY;
return Math.sqrt(dx * dx + dy * dy);
}
/**
* 数值计算梯度
*/
public double[] numericalGradient(double v0, double angle, double eps) {
double loss = lossFunction(v0, angle);
double lossVPlus = lossFunction(v0 + eps, angle);
double dLdv0 = (lossVPlus - loss) / eps;
double lossAPlus = lossFunction(v0, angle + eps);
double dLdAngle = (lossAPlus - loss) / eps;
return new double[]{dLdv0, dLdAngle};
}
/**
* 梯度下降寻找最优击球参数
* @param initV0 初始速度猜测
* @param initAngle 初始角度猜测（弧度）
* @param learningRate 学习率
* @param maxIterations 最大迭代次数
* @return 最优参数 [v0, angle]
*/
public double[] gradientDescent(double initV0, double initAngle,
double learningRate, int maxIterations) {
double v0 = initV0;
double angle = initAngle;
double eps = 0.001;
System.out.println("=== 梯度下降优化过程 ===");
System.out.printf("%-6s %-12s %-12s %-15s%n", "迭代", "速度", "角度(度)", "距离误差");
for (int i = 0; i < maxIterations; i++) {
double loss = lossFunction(v0, angle);
if (i % 50 == 0 || loss < holeRadius) {
System.out.printf("%-6d %-12.4f %-12.4f %-15.6f%n",
i, v0, Math.toDegrees(angle), loss);
}
if (loss < holeRadius) {
System.out.println("\n收敛！球落入洞中。");
break;
}
double[] grad = numericalGradient(v0, angle, eps);
// 参数更新
v0 -= learningRate * grad[0];
angle -= learningRate * grad[1];
// 约束：角度在 0~90 度之间
angle = Math.max(0.01, Math.min(Math.PI / 2 - 0.01, angle));
v0 = Math.max(1.0, v0); // 最小速度约束
}
return new double[]{v0, angle};
}
/**
* 二分搜索最优角度（固定速度）
* 当速度固定时，落点距离是角度的单峰函数
*/
public double findOptimalAngle(double v0, double targetX, double precision) {
double left = 0.01; // 最小角度（弧度）
double right = Math.PI / 2 - 0.01; // 最大角度
while (right - left > precision) {
double mid1 = left + (right - left) / 3;
double mid2 = right - (right - left) / 3;
double dist1 = Math.abs(simulateWithDrag(v0, mid1).x - targetX);
double dist2 = Math.abs(simulateWithDrag(v0, mid2).x - targetX);
if (dist1 < dist2) {
right = mid2;
} else {
left = mid1;
}
}
return (left + right) / 2;
}
/**
* 打印完整轨迹
*/
public void printTrajectory(double v0, double angle) {
List<double[]> trajectory = new ArrayList<>();
BallState state = new BallState(0, 0.01,
v0 * Math.cos(angle), v0 * Math.sin(angle));
while (state.y > GROUND_Y || state.vy > 0) {
trajectory.add(new double[]{state.x, state.y});
double v = Math.sqrt(state.vx * state.vx + state.vy * state.vy);
double dragAx = -AIR_RESISTANCE * v * state.vx;
double dragAy = -AIR_RESISTANCE * v * state.vy;
state.vx += dragAx * TIME_STEP;
state.vy += (dragAy - G) * TIME_STEP;
state.x += state.vx * TIME_STEP;
state.y += state.vy * TIME_STEP;
}
System.out.println("\n=== 轨迹采样 (每10步) ===");
System.out.printf("%-10s %-10s%n", "X(m)", "Y(m)");
for (int i = 0; i < trajectory.size(); i += 10) {
double[] p = trajectory.get(i);
System.out.printf("%-10.3f %-10.3f%n", p[0], p[1]);
}
BallState landing = simulateWithDrag(v0, angle);
System.out.printf("\n落点: (%.3f, %.3f) m%n", landing.x, landing.y);
System.out.printf("球洞: (%.3f, %.3f) m%n", holeX, holeY);
System.out.printf("误差: %.3f m%n", lossFunction(v0, angle));
}
/**
* 多初始点随机搜索 + 梯度下降
* 避免陷入局部最优
*/
public double[] randomRestartGradientDescent(int restarts, double learningRate, int iterations) {
Random rand = new Random(42);
double bestLoss = Double.MAX_VALUE;
double[] bestParams = null;
for (int r = 0; r < restarts; r++) {
double initV0 = 10 + rand.nextDouble() * 40; // 10~50 m/s
double initAngle = Math.toRadians(10 + rand.nextDouble() * 70); // 10~80度
double[] params = gradientDescent(initV0, initAngle, learningRate, iterations);
double loss = lossFunction(params[0], params[1]);
if (loss < bestLoss) {
bestLoss = loss;
bestParams = params;
}
}
return bestParams;
}
public static void main(String[] args) {
// 场景：球洞在100米外，半径0.1米
GolfSimulator golf = new GolfSimulator(100.0, 0.0, 0.1);
// 解析解（无阻力）
double v0Guess = Math.sqrt(G * 100 / Math.sin(2 * Math.toRadians(45)));
System.out.printf("解析猜测: v0=%.2f m/s, 角度=45度%n%n", v0Guess);
// 梯度下降优化
double[] optimal = golf.gradientDescent(v0Guess, Math.toRadians(45), 0.5, 500);
System.out.printf("\n最优参数: v0=%.4f m/s, 角度=%.4f度%n",
optimal[0], Math.toDegrees(optimal[1]));
// 打印轨迹
golf.printTrajectory(optimal[0], optimal[1]);
// 二分搜索验证
System.out.println("\n=== 二分搜索验证 ===");
double bisectAngle = golf.findOptimalAngle(optimal[0], 100.0, 0.0001);
System.out.printf("二分搜索最优角度: %.4f度%n", Math.toDegrees(bisectAngle));
}
}
算法复杂度分析
- 单次仿真：$O(T/\Delta t)$，$T$ 为飞行时间，$\Delta t$ 为时间步长。典型情况下约 $10^3$ 次迭代。
- 数值梯度计算：$O(T/\Delta t)$，需要对两个参数分别进行前向差分仿真。
- 梯度下降：$O(k \cdot T/\Delta t)$，$k$ 为迭代次数。
- 二分搜索：$O(\log(1/\epsilon) \cdot T/\Delta t)$，$\epsilon$ 为精度要求。
扩展思考
- 风速影响：可在物理模型中加入水平风力项 $F_w = c_w \cdot v_{wind}$，修正水平加速度。
- 地形起伏：将地面高度改为函数 $y = h(x)$，检测球与地面的碰撞点。
- 多杆策略：对于超长距离，可将问题建模为多阶段决策过程，每阶段选择使剩余距离最小的击球参数。
- 随机优化：实际球场存在地面硬度、草皮摩擦等不确定性，可引入高斯噪声进行鲁棒优化。
总结
高尔夫击球问题展示了如何将连续优化与物理仿真相结合。梯度下降通过数值差分逼近目标函数的局部最优点，而多初始点随机重启策略有效降低了陷入局部最优的风险。抛物线物理模型作为基础框架，可通过增加阻力、风力等修正项不断逼近真实场景。
