CODEWITHSUNDEEP
pythonoptimizationodegekko
Aleksejs Fomins
Aleksejs Fomins

Reputation: 900

How to fix Gekko's `@error: Solution not found`

I want to perform parameter inference for an ODE model (m.options.IMODE = 5). I have 2 models in Gekko. First, I generate model data from each model, and perform parameter inference directly on the corresponding model data. That works no problem for each model. Then, I attempt to perform parameter inference on real data, that has roughly the same number of datapoints as the model data. The first model finds a reasonable fit, the second one crashes with @error: Solution not found. I would like to learn the procedure of debugging problems like this.

Firstly, the last few lines of the Gekko error report for m.solve(disp=True) is

 Iter    Objective  Convergence
  250  8.46435E+03  2.23511E-01
 Maximum iterations
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :    101.300599999988      sec
 Objective      :    8182.67185576400     
 Unsuccessful with error code            0
 ---------------------------------------------------
 
 Creating file: infeasibilities.txt
 Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
 @error: Solution Not Found

Question 1: Is it correct to assume that the solver terminated because it reached the maximum number of iterations, and not because if infeasibilities?

I have tried all of the suggestions from this answer regarding reaching maximum iteration number, but to no avail.

  1. I have tried random and fixed (somewhat reasonable) initialization for unknown parameters. In all cases, the objective function seems stuck at the exact same value on every iteration. I presume that increasing the number of iterations would not fix this problem.
  2. I have tried APOPT, BPOPT and IPOPT. For APOPT and IPOPT, the solver behaves the same, crashing with Solution Not Found having made no progress on the objective function. With BPOPT, it crashes with WARNING MESSAGE FROM SUBROUTINE MA27BD *** INFO(1) = 3 MATRIX IS SINGULAR.
  3. I have tried using cold start as described in the mentioned answer. For all solvers, I get Successful Pre-solve Solution for the cold start, and still the Solution Not Found for the second part.

Question 2: How do I proceed to interpret what is going wrong and fix this problem?

I'm not allowed to share the original data, but I can share the ODE. For the first model (which fits fine):

m.Equation(A * x.dt() == K * B * y * (x1 - x) + C * (x2 - x))

z_est = x + (x1 - x) * (1 - B)
objective = m.Intermediate((z - z_est)**2)
m.Minimize(objective)

Here A, B, C are the unknown parameters we seek to infer, K is a known constant, and y, x1, x2, z are known time-dependent continuous variables, which are measured for every timestep. The only difference to the second model (that does not fit is)

B = m.Intermediate(m.exp(-y0 / y))

Now, the unknown parameters are A, y0, C. There are no additional inequality constraints. The only constraints are specifying very broad upper and lower bounds on the unknown parameters, such as

bounds = dict(A = [0.0001, 10000000], C = [0.0001, 100000], y0 = [0.001, 1000])

Question 3: Is there any use in looking at the infeasibilities file? I have looked at it briefly for APOPT, and I have found that POSSIBLE INFEASBILE EQUATIONS section is empty, but there are a lot of complicated lines in the ACTIVE OBJECTIVE EQUATIONS, ACTIVE EQUATIONS and INACTIVE EQUATIONS. Some of the lines in ACTIVE EQUATIONS say not available at the end.

Upvotes: 2

Views: 79

Answers (1)

John Hedengren
John Hedengren

Reputation: 14346

The problem may be with divide-by-zero as the solver is searching for a solution. Here is the equation that may be problematic:

B = m.Intermediate(m.exp(-y0 / y))

A few strategies to deal with the problem are:

  1. Rearrange equation to avoid divide-by-zero
B = m.Var(lb=1e-3)
m.Equation(y*m.log(B)===-y0)
  1. Add a constraint to y
y = m.Var(lb=1e-3) # at initialization
y.LOWER = 1e-3     # or set bound independently
  1. Add a small quantity to the denominator
y = m.Var(lb=0)
B = m.Intermediate(m.exp(-y0 / (y+1e-3)))

However, the third option changes the equation slightly.

If none of these suggestions help, I recommend creating a fixed seed random data set that reproduces the issue. It is much easier to debug a problem when there is a minimal example that demonstrates the issue. If it is an initialization problem, there are additional suggestions in this paper:

Safdarnejad, S.M., Hedengren, J.D., Lewis, N.R., Haseltine, E., Initialization Strategies for Optimization of Dynamic Systems, Computers and Chemical Engineering, 2015, Vol. 78, pp. 39-50, DOI: 10.1016/j.compchemeng.2015.04.016.

Upvotes: 0

Related Questions