GEKKO NL Solver拟合值与参考值不匹配问题排查：双参数模型拟合失效，单参数模型表现更优

Let's break down why your two-parameter model is spitting out near-zero predictions and how to fix this issue:

Possible Root Causes #

1. Stuck in a Local Optimum #

Your two-parameter model has more degrees of freedom, which makes it easier for the solver to land in a unhelpful local minimum. In this case, the solver might be choosing a1 and a2 values that make the denominator 1 - (a1*x)² + (a2*x)⁴ extremely large—this drives y close to zero, even though it's a terrible fit for your measured data. The single-parameter model is simpler, so it's far less likely to get trapped in such a bad local optimum.

2. Unconstrained Parameter Ranges #

You didn't set any bounds on a1 and a2, so the solver can pick values that blow up the denominator. For example, a large a2 will make (a2*x)^4 dominate the denominator, turning it into a huge number and y into a tiny one.

3. Suboptimal Initial Guess #

While your 0.3 initial guess prevents solver errors, it might not be close enough to the true parameter values to guide the solver toward a good global optimum.

Fixes to Try #

1. Add Parameter Bounds #

Restrict a1 and a2 to reasonable ranges based on your data. For example, if your x values are small, you could set bounds like 0 <= a1 <= 1 and 0 <= a2 <= 0.5 to keep the denominator from exploding.

2. Reparameterize the Model #

Rewrite the model to use non-negative parameters for the polynomial terms—this makes optimization more stable and aligns with physical intuition (if your data comes from a real-world system):

# Use b1 = a1² and b2 = a2⁴ (both non-negative)
b1 = m.FV(value=0.09, lb=0)  # 0.3² = 0.09
b2 = m.FV(value=0.0081, lb=0) # 0.3⁴ = 0.0081
b1.STATUS = 1
b2.STATUS = 1

# Update the equation to use the new parameters
m.Equation(y == x / (1 - b1*x**2 + b2*x**4))

3. Enable Global Optimization #

GEKKO's APOPT solver supports global search to avoid local minima. Enable it with these settings:

m.options.SOLVER = 1  # Use APOPT solver
m.options.MAX_ITER = 1000
m.solver_options = ['minlp_maximum_iterations 1000',
                    'minlp_gap_tol 1.0e-4',
                    'minlp_integer_tol 1.0e-4']

4. Normalize Your Data #

If your x or y values are very large or small, normalize them to the [0,1] range to improve numerical stability:

# Normalize data
xm_norm = (xm - xm.min()) / (xm.max() - xm.min())
ym_norm = (ym - ym.min()) / (ym.max() - ym.min())

# Use normalized values in the model
x = m.Param(value=xm_norm)
y = m.CV(value=ym_norm)
y.FSTATUS = 1

After optimization, reverse the normalization to get predictions in the original scale.

Improved Code Example #

Here's a revised version of your code with parameter bounds, reparameterization, and global optimization enabled:

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

# Load measured data
measure_data = np.load("C:/measure.npy")
xm = measure_data[0]
ym = measure_data[1]

# GEKKO model setup
m = GEKKO()

# Normalize data for better numerical stability
xm_norm = (xm - xm.min()) / (xm.max() - xm.min())
ym_norm = (ym - ym.min()) / (ym.max() - ym.min())

# Reparameterized parameters (non-negative)
x = m.Param(value=xm_norm)
b1 = m.FV(value=0.09, lb=0, ub=2)  # b1 = a1²
b2 = m.FV(value=0.0081, lb=0, ub=1) # b2 = a2⁴
b1.STATUS = 1
b2.STATUS = 1

# Variables
y = m.CV(value=ym_norm)
y.FSTATUS = 1

# Model equation
m.Equation(y == x / (1 - b1*x**2 + b2*x**4))

# Regression mode with global optimization
m.options.IMODE = 2
m.options.SOLVER = 1  # APOPT solver
m.options.MAX_ITER = 1000
m.solver_options = ['minlp_maximum_iterations 1000',
                    'minlp_gap_tol 1.0e-4']

# Run optimization (enable disp to see solver output)
m.solve(disp=True)

# Convert back to original a1, a2 parameters
b1_opt = b1.value[0]
b2_opt = b2.value[0]
a1_opt = np.sqrt(b1_opt)
a2_opt = np.power(b2_opt, 1/4)

# Generate predictions and reverse normalization
optimized_y_norm = xm_norm / (1 - b1_opt*xm_norm**2 + b2_opt*xm_norm**4)
optimized_y = optimized_y_norm * (ym.max() - ym.min()) + ym.min()

# Plot results (include single-parameter fit for comparison)
plt.figure(1)
plt.plot(xm, ym, 'k', label="Measurements")
plt.plot(xm, optimized_y, 'r', label="Optimized (Two-Parameter)")

# Add single-parameter model fit
m_single = GEKKO()
x_single = m_single.Param(value=xm_norm)
y_single = m_single.CV(value=ym_norm)
y_single.FSTATUS = 1
a1_single = m_single.FV(value=0.3)
a1_single.STATUS = 1
m_single.Equation(y_single == x_single / (1 - (a1_single*x_single)**2))
m_single.options.IMODE = 2
m_single.solve(disp=False)

optimized_y_single_norm = xm_norm / (1 - (a1_single.value[0]*xm_norm)**2)
optimized_y_single = optimized_y_single_norm * (ym.max() - ym.min()) + ym.min()
plt.plot(xm, optimized_y_single, 'b--', label="Optimized (Single-Parameter)")

plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

print(f"Optimized a1: {a1_opt:.4f}, a2: {a2_opt:.4f}")

Key Notes #

• Check the solver output (enable disp=True) to confirm it converged to a feasible solution.
• If the two-parameter model still doesn't outperform the single-parameter one, it might mean your data doesn't actually need the quartic term—sometimes simpler models are the best fit!

内容的提问来源于stack exchange，提问作者audi02