Question: I want to do parameter estimation on a Spring-Mass system  with direct collocation method. The parameter k should be determined from response. I simulated this system by from scipy.integrate import odeint import numpy as np def dy(y, t): x, xdot = y return [xdot, -50*x] t = np.linspace(0, 1, 40) sol = odeint(dy, [2.0, 1.0], t) sol_x = sol[:, 0] sol_xdot = sol[:, 1]  Then I have these code to identify parameter k: from gekko import GEKKO m = GEKKO(remote=False) m.time = t x = m.CV(value=sol_x); x.FSTATUS = 1 # fit to measurement xdot = m.CV(value=sol_xdot); xdot.FSTATUS = 1 k = m.FV(value = 40.0); k.STATUS = 1 # change initial value of k here m.Equation(x.dt() == xdot) # differential equation m.Equation(xdot.dt() == -k*x) m.options.IMODE = 5 # dynamic estimation m.options.NODES = 40 # collocation nodes m.options.EV_TYPE = 2 m.solve(disp=False) # display solver output  By playing around initial value of k, I found k will converge to real value 50 if its initial value is 25 to 65. Otherwise the result will be -0.39 which is not good. I'm quite confusing because this system is linear and should be easy to be solved. So my quietion: how to fine tune the above code so that k converge to 50 with arbitry initial value?

Answer: The -0.39 is a local minimum to the optimization problem. As the initial guess is further from the solution, it finds a different local solution. To prevent non-physical solutions, add a constraint for the solver to search only within bounds. This can be done during initialization with: k = m.FV(value=ki,lb=10,ub=100)  or after initialization with: k.LOWER = 10 k.UPPER = 100  Here is a complete script that returns the correct solution regardless of the initial guess. from gekko import GEKKO from scipy.integrate import odeint import numpy as np # generate data def dy(y, t): x, xdot = y return [xdot, -50*x] t = np.linspace(0, 1, 40) sol = odeint(dy, [2.0, 1.0], t) sol_x = sol[:, 0] sol_xdot = sol[:, 1] # regression with range of initial guess values kx = np.linspace(0,100,11) for i,ki in enumerate(kx): m = GEKKO(remote=False) m.time = t x = m.CV(value=sol_x); x.FSTATUS = 1 xdot = m.CV(value=sol_xdot); xdot.FSTATUS = 1 k = m.FV(value=ki,lb=10,ub=100); k.STATUS = 1 m.Equation(x.dt() == xdot) m.Equation(xdot.dt() == -k*x) m.options.IMODE = 5 # dynamic estimation m.options.NODES = 6 # collocation nodes (2-6) m.options.EV_TYPE = 2 m.solve(disp=False) print(f'Initial Guess: {ki} ' + f'Solution: {k.value[0]} ' + f'ObjFcn: {m.options.OBJFCNVAL}')  The output shows that the correct solution is found, regardless of the initial guess. Initial Guess: 0.0 Solution: 50.000000284 ObjFcn: 1.8231507179e-10 Initial Guess: 10.0 Solution: 50.000000284 ObjFcn: 1.8229910145e-10 Initial Guess: 20.0 Solution: 50.000000284 ObjFcn: 1.8231237768e-10 Initial Guess: 30.0 Solution: 50.000000284 ObjFcn: 1.8229796958e-10 Initial Guess: 40.0 Solution: 50.000000284 ObjFcn: 1.8229906845e-10 Initial Guess: 50.0 Solution: 50.000000284 ObjFcn: 1.8229906842e-10 Initial Guess: 60.0 Solution: 50.000000284 ObjFcn: 1.8229906933e-10 Initial Guess: 70.0 Solution: 50.000000284 ObjFcn: 1.8229906849e-10 Initial Guess: 80.0 Solution: 50.000000284 ObjFcn: 1.8230214413e-10 Initial Guess: 90.0 Solution: 50.000000284 ObjFcn: 1.8229908716e-10 Initial Guess: 100.0 Solution: 50.000000284 ObjFcn: 1.8230004454e-10  Without the constraints, a local solution is found but the objective function is higher that indicates it is not a good fit. Initial Guess: 0.0 Solution: -0.39654383084 ObjFcn: 88263.987254 Initial Guess: 10.0 Solution: -0.39654383084 ObjFcn: 88263.987254 Initial Guess: 20.0 Solution: -0.39654383084 ObjFcn: 88263.987254 Initial Guess: 30.0 Solution: 50.000000284 ObjFcn: 1.8229906872e-10 Initial Guess: 40.0 Solution: 50.000000284 ObjFcn: 1.8229906866e-10 Initial Guess: 50.0 Solution: 50.000000284 ObjFcn: 1.8229906856e-10 Initial Guess: 60.0 Solution: 50.000000284 ObjFcn: 1.8229906861e-10 Initial Guess: 70.0 Solution: -0.39654383084 ObjFcn: 88263.987254 Initial Guess: 80.0 Solution: -0.39654383084 ObjFcn: 88263.987254 Initial Guess: 90.0 Solution: -0.39654383084 ObjFcn: 88263.987254 Initial Guess: 100.0 Solution: -0.39654383085 ObjFcn: 88263.987254  If constraints are unknown or there are multiple local solutions then a multi-start method can be used to search for a global optimum as shown in the Design Optimization course on Global Optimization. Below is a Global search over an initial guess space. The hyperopt package uses Bayesian optimization to find the global solution. from gekko import GEKKO from scipy.integrate import odeint import numpy as np from hyperopt import fmin, tpe, hp from hyperopt import STATUS_OK, STATUS_FAIL # generate data def dy(y, t): x, xdot = y return [xdot, -50*x] t = np.linspace(0, 1, 40) sol = odeint(dy, [2.0, 1.0], t) sol_x = sol[:, 0] sol_xdot = sol[:, 1] def objective(params): ki = params['kx'] m = GEKKO(remote=False) m.time = t x = m.CV(value=sol_x); x.FSTATUS = 1 xdot = m.CV(value=sol_xdot); xdot.FSTATUS = 1 k = m.FV(value=ki); k.STATUS = 1 m.Equation(x.dt() == xdot) m.Equation(xdot.dt() == -k*x) m.options.IMODE = 5 # dynamic estimation m.options.NODES = 6 # collocation nodes (2-6) m.options.EV_TYPE = 2 m.solve(disp=False) obj = m.options.OBJFCNVAL if m.options.APPSTATUS==1: s=STATUS_OK else: s=STATUS_FAIL m.cleanup() return {'loss':obj, 'status': s, 'k':k.value[0]} # Define the search space for the hyperparameters space = {'kx': hp.quniform('kx', 0, 100, 10),} best = fmin(objective, space, algo=tpe.suggest, max_evals=10) sol = objective(best) print(f"Solution Status: {sol['status']}") print(f"Objective: {sol['loss']:.2f}") print(f"Initial Guess: {best['kx']}") print(f"Solution: {sol['k']}")  Here is the solution: Solution Status: ok Objective: 0.00 Initial Guess: 30.0 Solution: 50.000000284  While not needed for this problem, there is logic to detect when the poor initial guess produces a failed solution and eliminates that initial guess.