here comes a chemical problem statement, describing the kinetics of an reacting mixture between CO2 and H2 with increasing T. I would be glad to get some hints on how to improve this code to solve a NLP. As depicted, I still face non-physical results, in this case negative concentration values for CO mainly. Any suggestions, please? Thank you. And if you decide to let run this with statements like: m.Equation(0 == -(k1 * S * CO2 * H2 * (1 - ((CO * H2O) / (CO2 * H2)) / K1eq_param))) it will give a much better fit, but with only zero-noisy values for the kinetic rate konstants k import numpy as np from gekko import GEKKO import matplotlib.pyplot as plt import csv CH4_meas = np.array([0.1361,0.1885,0.2481,0.3341,0.4266]) # and further values... CO_meas = np.array([0.0281,0.0376,0.0476,0.0628,0.0817]) H2O_meas = np.array([0.1201,0.234,0.3821,0.5821,0.8274]) H2_meas = np.array([62.6176,62.3346,62.0223,61.7904,61.5901]) CO2_meas = np.array([30.4669,30.3408,30.3541,30.3541,30.2984]) K1eq = np.array([0.01690,0.01881,0.02088,0.02312,0.02555]) K2eq = np.array([5.7990E+08,3.3957E+08,2.0132E+08,1.2079E+08,7.3315E+07]) K3eq = np.array([9.7992E+06,6.3860E+06,4.2031E+06,2.7930E+06,1.8732E+06]) CH4_values = [] CO_values = [] H2O_values = [] H2_values = [] CO2_values = [] k1_values = [] k2_values = [] k3_values = [] S=25 for i in range(len(CH4_meas)): m = GEKKO(remote=False) # Initial guess for the rate constants, but when I set lb=0, no solution can be found k1 = m.FV(value=0.0025, lb=None, ub=None) k2 = m.FV(value=0.0025, lb=None, ub=None) k3 = m.FV(value=0.2, lb=None, ub=None) k1.STATUS = 1 k2.STATUS = 1 k3.STATUS = 1 # Controlled Variables (CV) with measurement normalization CH4 = m.CV(CH4_meas / (H2_meas[0] + CO2_meas[0])) CO = m.CV(CO_meas / (H2_meas[0] + CO2_meas[0])) H2O = m.CV(H2O_meas / (H2_meas[0] + CO2_meas[0])) H2 = m.CV(H2_meas / (H2_meas[0] + CO2_meas[0])) CO2 = m.CV(CO2_meas / (H2_meas[0] + CO2_meas[0])) CH4.FSTATUS = 1 CO.FSTATUS = 1 H2O.FSTATUS = 1 H2.FSTATUS = 1 CO2.FSTATUS = 1 # Parameters K1eq_param = m.Param(value=K1eq) K2eq_param = m.Param(value=K2eq) K3eq_param = m.Param(value=K3eq) # Steady-state equations m.Equation(CO2 == -(k1 * S * CO2 * H2 * (1 - ((CO * H2O) / (CO2 * H2)) / K1eq_param)) - k3 * S * CO2 * (1 - ((CH4 * H2O**2) / (CO2 * H2**4)) / K3eq_param)) m.Equation(H2 == -(k1 * S * CO2 * H2 * (1 - ((CO * H2O) / (CO2 * H2)) / K1eq_param) - 3 * k2 * S * CO * (1 - ((CH4 * H2O) / (CO * H2**3)) / K2eq_param) - 4 * k3 * S * CO2 * (1 - ((CH4 * H2O**2) / (CO2 * H2**4)) / K3eq_param))) m.Equation(CO == (k1 * S * CO2 * H2 * (1 - ((CO * H2O) / (CO2 * H2)) / K1eq_param) - k2 * S * CO * (1 - ((CH4 * H2O) / (CO * H2**3)) / K2eq_param))) m.Equation(H2O == (k1 * S * CO2 * H2 * (1 - ((CO * H2O) / (CO2 * H2)) / K1eq_param) + k2 * S * CO * (1 - ((CH4 * H2O) / (CO * H2**3)) / K2eq_param) + k3 * S * CO2 * (1 - ((CH4 * H2O**2) / (CO2 * H2**4)) / K3eq_param))) m.Equation(CH4 == (k2 * S * CO * (1 - ((CH4 * H2O) / (CO * H2**3)) / K2eq_param) + k3 * S * CO2 * (1 - ((CH4 * H2O**2) / (CO2 * H2**4)) / K3eq_param))) # Objective function for regression fitting m.Obj((CH4.value - CH4_meas)**2 + (CO.value - CO_meas)**2 + (H2O.value - H2O_meas)**2 + (H2.value - H2_meas)**2 + (CO2.value - CO2_meas)**2) # Set the mode for steady-state optimization m.options.IMODE = 3 # Steady-state optimization m.options.EV_TYPE = 1 #2 = Sum of squared error, 1 = Sum of absoluted error m.options.MAX_ITER = 300 m.options.OTOL = 1e-6 m.options.RTOL = 1e-6 # Solve the model m.solve(disp=True) # Store results CH4_values.append(CH4.value[0]) CO_values.append(CO.value[0]) H2O_values.append(H2O.value[0]) H2_values.append(H2.value[0]) CO2_values.append(CO2.value[0]) k1_values.append(k1.value[0]) k2_values.append(k2.value[0]) k3_values.append(k3.value[0]) print("k1 = " + str(k1.value[0])) print("k2 = " + str(k2.value[0])) print("k3 = " + str(k3.value[0])) with open('OptimGekko_Methane_RWGS.csv', 'w', newline='') as csvfile: csvwriter = csv.writer(csvfile) csvwriter.writerow(['T', 'k1', 'k2', 'k3','CO2','H2','CO','H2O','CH4']) for i in range(len(k1_values)): csvwriter.writerow([time, k1_values, k2_values, k3_values, CO2_values, H2_values,CO_values, H2O_values[i], CH4_values[i]]) plt.figure(figsize=(10, 12)) plt.subplot(5, 1, 1) plt.plot(time[:len(CH4_values)], CH4_values, label='CH4_simul') plt.plot(time[:len(CH4_meas)], CH4_meas/(30.47+62.62), 'ro', label='CH4_meas') plt.ylabel('c/c0 (-)') plt.xlabel('T (°C)') plt.legend() plt.grid() plt.title('CH4') plt.subplot(5, 1, 2) plt.plot(time[:len(CO_values)], CO_values, label='CO_simul') plt.plot(time[:len(CO_meas)], CO_meas/(30.47+62.62), 'ro', label='CO_meas') plt.ylabel('c/c0 (-)') plt.xlabel('T (°C)') plt.legend() plt.grid() plt.title('CO') plt.subplot(5, 1, 3) plt.plot(time[:len(H2O_values)], H2O_values, label='H2O_simul') plt.plot(time[:len(H2O_meas)], H2O_meas/(30.47+62.62), 'ro', label='H2O_meas') plt.ylabel('c/c0 (-)') plt.xlabel('T (°C)') plt.legend() plt.grid() plt.title('H2O') plt.subplot(5, 1, 4) plt.plot(time[:len(H2_values)], H2_values, label='H2_simul') plt.plot(time[:len(H2_meas)], H2_meas/(30.47+62.62), 'ro', label='H2_meas') plt.ylabel('c/c0 (-)') plt.xlabel('T (°C)') plt.legend() plt.grid() plt.title('H2') plt.subplot(5, 1, 5) plt.plot(time[:len(CO2_values)], CO2_values, label='CO2_simul') plt.plot(time[:len(CO2_meas)], CO2_meas/(30.47+62.62), 'ro', label='CO2_meas') plt.ylabel('c/c0 (-)') plt.xlabel('T (°C)') plt.legend() plt.grid() plt.title('CO2') plt.tight_layout() plt.savefig('Optim_RWGS_CH4.png') plt.show() [url="https://stackoverflow.com/questions/79271590/how-to-avoid-non-physical-results-in-solving-a-nlp"]Continue reading...[/url][/i][/i]
