"Exception: @error: Solution not Found" on and off while solving the same problem

I am solving multiple optimization problems with the Gekko optimization suite. I generate random numbers for three variables (book3.csv) and import them to python as initial values and parameters. I run optimization solver and collect the results. It works well. However, when I tried this about 20 times to check, I received an error message of Exception: @error: Solution not Found several times. It stops at the same location: ...gekko.py line 2174, in solve raise Exception(response)

I collected the value of control variables, u, at each round of optimization, i, at each time points, j, and saved it at a dataframe, utot, in order to track all the decisions made at each time for each round. I noticed that the last column (at the last time points) of utot has some positive values, which I think shouldn't be. I intentionally put Pe0=0, so that u=0 is optimum at all time points, but still the last column of utot has some positive value. I am not sure what am I missing here.

from gekko import GEKKO
import pandas as pd
import numpy as np
from numpy import random


file_path='./Book3.csv'
d=pd.read_csv("Book3.csv")
dftemp = pd.DataFrame(data=d)
na=len(dftemp)


# time points
n=51 
year=50

# constants
Pa0 = 1 
Pe0 = 0
C = 10 
r = 0.05 
k=50
shift=100
ll=0.5*random.rand(n)

index = range(0, na)
columns = range(0, n)
columns1 = range(0, 1)

#Dataframe to collect results
Atot=pd.DataFrame(index=index,columns=columns)
Etot=pd.DataFrame(index=index,columns=columns)
utot=pd.DataFrame(index=index,columns=columns)
PVtot=pd.DataFrame(index=index,columns=columns1)
    
for i in range(0,na):
    # create GEKKO model
    m = GEKKO(remote=False)

    m.time = np.linspace(0,year,n)
    t=m.time

    
    A0=dftemp.loc[i][0]
    h=dftemp.loc[i][1]
    emax=dftemp.loc[i][2]    

    u = m.MV(value=0,lb=0, ub=emax) 
    u.STATUS = 1
    u.DCOST = 0
    
    A = m.SV(value=A0)
    E = m.SV(value=0)
    t = m.Param(value=m.time)
    Pe = m.Var(value=Pe0)
    d = m.Var(value=1)
    l = m.Param(value=ll)
    
    # Equation
    m.Equation(A.dt()==-u)
    m.Equation(E.dt()==u)
    m.Equation(Pe.dt()==-Pe/k)
    m.Equation(d==m.exp(-t*r))
    m.Equation(A>=0)
        
    # Objective (Utility)
    J = m.Var(value=0)
    
    # Final objective
    Jf = m.FV()
    Jf.STATUS = 1
    m.Connection(Jf,J,pos2='end')
    m.Equation(J.dt() == m.log(A+E*(1-l))*h*Pa0-r*C*u+E*Pe+shift)*d
    
  
    # maximize U
    m.Maximize(Jf)
    
    # options
    m.options.IMODE = 6  # optimal control
    m.options.NODES = 3  # collocation nodes
    m.options.SOLVER = 3 # solver (IPOPT)
    
    # solve optimization problem
    m.solve()
    
    # print profit
    print('Optimal Profit: ' + str(Jf.value[0]))
 

    for j in range(0,n):
        Atot.loc[i][j]=A.value[j]
        Etot.loc[i][j]=E.value[j]
        utot.loc[i][j]=u.value[j]
    PVtot.loc[i][0]=Jf.value[0]
    
print(PVtot.sum())