问题

I'm currently implementing the Galerkin-Method for finite Elements by hand and using python gekko to solve the resulting non-linear algebraic equation system. This creates no issues for small systems and works fine.

As soon as the system becomes more complex involving long equations and exponential terms with m.exp() the equations might not become formatted properly anymore for the solver. The exact error message is @error: Equation Definition Equation without an equality (=) or inequality (>,<).

I checked the equations via m.open_folder() in the model.apm. And it seems, that somehow a line break is added into the equation loosing the mathematical relationship between the resulting expressions. This is why the solver does not recognize equality (=) or inequality (>,<) in the first expression.
Also I found out that the exponential term m.exp(-EA / (R * T)) causes the line break. Replacing it with np.exp(-EA/(R*T0)) hotfixes the issue. It might be due to the fact the gekko shortens long equations in the background and there is a bug within this procedure. Or something is wrong with my equation definition. Here is my code.

Basis for generating long and complex algebraic equations

import pandas as pd
import numpy as np
from gekko import GEKKO
from typing import List, Callable
m = GEKKO(remote=False)
# constants
# -------------------------
R = 8.314  # J/mol K
# Rahmenbedingungen
# -------------------------
T0 = 550  # K
p0 = 15 * 10 ** 5 # Pa
# components
# -------------------------
x_bulk = np.array(
    [0.1, 0.08, 0.05]
)  # no unit Ethen, O2, CO2, EO, H2O
# catalyst pellet
# -------------------------
rho_p = 1171.2  # kg / m**3 Dichte eines Pellets
epsilon_p = 0.7  # Poroesitaet
tau = 4  # Tortuositaet
# heat transfer
# -------------------------
lambda_pm = 16  # W/m K Wärmeleitfähigkeit des reinen Katalysatorfestoffes
lambda_p = (1 - epsilon_p) * lambda_pm  # W/m K Wärmeleitfähigkeit des porösen Pellets
# reaction enthalpies
# -------------------------
dh_r1 = -107000  # J/mol
dh_r2 = -1323000  # J/mol
Deff = 5.42773893 * 10 **7

# reaction speeds 
# --------
def r_r1(p: np.array, T: float = T0) -> float:
    n_E = 0.6  # Reaktionsordnung Ethen
    n_O2 = 0.5  # Reaktionsordnung Sauerstoff
    k0 = 6.275 * 10**6  # pre-exponential collision factor [mol / (kg s Pa ** (nE+nO2))]
    EA = 74900 # activation energy [J/mol]
    K0 = 1.985 * 10**2  # pre-exponential adsorption factor [1/Pa]
    Tads = 2400  # adsorption temperature [K]
    r = k0 * m.exp(-EA / (R * T)) * p[0] ** n_E * p[1] ** n_O2 / (1 + K0 * np.exp(Tads / T0) * p[2])
    # this is m.exp() one cause for the quations not being formatted correctly
    # change it to np.exp(-EA/(R*T0)) and the equations are formatted correctly for the solver
    # this may be a mistake in my code or a bug in gekko
    return r

# differential equations 
# --------
def rhs(x: float, y: np.array, dy: np.array) -> np.array:
    # second derivatives of the partial pressures [bar/m**2]
    # use bar instead of Pa to increase numerical precision
    rhs_E   = -2 / x * dy[0] + R * y[-1] * T0 * rho_p / (Deff * p0) * (r_r1(p=y[0:-1] * p0, T=y[-1] * T0))
    rhs_O2  = -2 / x * dy[1] + R * y[-1] * T0 * rho_p / (Deff * p0) * (0.5 * r_r1(p=y[0:-1] * p0, T=y[-1] * T0))
    rhs_CO2 = -2 / x * dy[2] + R * y[-1] * T0 * rho_p / (Deff * p0) * (-2) * r_r1(p=y[0:-1] * p0, T=y[-1] * T0)
    rhs_T   = -2 / x * dy[-1] - rho_p / (lambda_p * T0) * (r_r1(p=y[0:-1] * p0, T=y[-1] * T0) * (-dh_r1))
    return np.array([rhs_E, rhs_O2, rhs_CO2, rhs_T])

# quadrature over the right hand side of the difffereantial equation
# ----------------

def give_quad_expressions(
    rhs: Callable, phi: List[np.poly1d], dphi: List[np.poly1d]
) -> Callable:
    """Creates an expression for the quadrature of the right-hand-side of the
    differential equation that is only dependent on the function values at the
    nodes y_i. It takes the quadrature base points to evaluate the
    base functions phi_i and their derivatives dphi_i to set up the expression.

    Args:
        rhs (Callable): Right hand side of the differential equation.
        Needs to be a function of x, y, and dy
        phi (List[np.poly1d]): List with all base polynomials phi. Mostly, base polynomials of Legendre type.
        dphi (List[np.poly1d]): List with the corresponding derivatives of the base polynomials.

    Returns:
        Callable: Expression for the weighted quadratures of the right-hand-side for each node that is only
        dependent on the unknown function values y_i. It can be called by a solving algorithm.
    """
    # these will be defined during runtime
    number_base_points = 2
    base_points = np.array([0.00021132, 0.00078868])
    weights = np.array([0.0005, 0.0005])

    # evaluate base functions at quadrature base points
    phi_eval = np.zeros(shape=(len(phi), number_base_points))
    dphi_eval = np.zeros(shape=(len(dphi), number_base_points))
    for index, phi_i in enumerate(phi):
        phi_eval[index] = phi_i(base_points)
        dphi_eval[index] = dphi[index](base_points)

    # set

<details>
<summary>英文:</summary>

I&#39;m currently implementing the Galerkin-Method for finite Elements by hand and using python gekko to solve the resulting non-linear algebraic equation system. This creates no issues for small systems and works fine.

As soon as the system becomes more complex involving long equations and exponential terms with `m.exp()` the equations might not become formatted properly anymore for the solver. The exact error message is `@error: Equation Definition Equation without an equality (=) or inequality (&gt;,&lt;)`.

I checked the equations via `m.open_folder()` in the `model.apm`. And it seems, that somehow a line break is added into the equation loosing the mathematical relationship between the resulting expressions. This is why the solver does not recognize `equality (=) or inequality (&gt;,&lt;)` in the first expression.
Also I found out that the exponential term `m.exp(-EA / (R * T))` causes the line break. Replacing it with `np.exp(-EA/(R*T0))` hotfixes the issue. It might be due to the fact the gekko shortens long equations in the background and there is a bug within this procedure. Or something is wrong with my equation definition. Here is my code.

Basis for generating long and complex algebraic equations

import pandas as pd
import numpy as np
from gekko import GEKKO
from typing import List, Callable
m = GEKKO(remote=False)

constants

-------------------------

R = 8.314 # J/mol K

Rahmenbedingungen

-------------------------

T0 = 550 # K
p0 = 15 * 10 ** 5 # Pa

components

-------------------------

x_bulk = np.array(
[0.1, 0.08, 0.05]
) # no unit Ethen, O2, CO2, EO, H2O

catalyst pellet

-------------------------

rho_p = 1171.2 # kg / m**3 Dichte eines Pellets
epsilon_p = 0.7 # Poroesitaet
tau = 4 # Tortuositaet

heat transfer

-------------------------

lambda_pm = 16 # W/m K Wärmeleitfähigkeit des reinen Katalysatorfestoffes
lambda_p = (1 - epsilon_p) * lambda_pm # W/m K Wärmeleitfähigkeit des porösen Pellets

reaction enthalpies

-------------------------

dh_r1 = -107000 # J/mol
dh_r2 = -1323000 # J/mol
Deff = 5.42773893 * 10 **7

reaction speeds

--------

def r_r1(p: np.array, T: float = T0) -> float:
n_E = 0.6 # Reaktionsordnung Ethen
n_O2 = 0.5 # Reaktionsordnung Sauerstoff
k0 = 6.275 * 106 # pre-exponential collision factor [mol / (kg s Pa ** (nE+nO2))]
EA = 74900 # activation energy [J/mol]
K0 = 1.985 * 102 # pre-exponential adsorption factor [1/Pa]
Tads = 2400 # adsorption temperature [K]
r = k0 * m.exp(-EA / (R * T)) * p[0] ** n_E * p[1] ** n_O2 / (1 + K0 * np.exp(Tads / T0) * p[2])
# this is m.exp() one cause for the quations not being formatted correctly
# change it to np.exp(-EA/(R*T0)) and the equations are formatted correctly for the solver
# this may be a mistake in my code or a bug in gekko
return r

differential equations

--------

def rhs(x: float, y: np.array, dy: np.array) -> np.array:
# second derivatives of the partial pressures [bar/m**2]
# use bar instead of Pa to increase numerical precision
rhs_E = -2 / x * dy[0] + R * y[-1] * T0 * rho_p / (Deff * p0) * (r_r1(p=y[0:-1] * p0, T=y[-1] * T0))
rhs_O2 = -2 / x * dy[1] + R * y[-1] * T0 * rho_p / (Deff * p0) * (0.5 * r_r1(p=y[0:-1] * p0, T=y[-1] * T0))
rhs_CO2 = -2 / x * dy[2] + R * y[-1] * T0 * rho_p / (Deff * p0) * (-2) * r_r1(p=y[0:-1] * p0, T=y[-1] * T0)
rhs_T = -2 / x * dy[-1] - rho_p / (lambda_p * T0) * (r_r1(p=y[0:-1] * p0, T=y[-1] * T0) * (-dh_r1))
return np.array([rhs_E, rhs_O2, rhs_CO2, rhs_T])

quadrature over the right hand side of the difffereantial equation

----------------

def give_quad_expressions(
rhs: Callable, phi: List[np.poly1d], dphi: List[np.poly1d]
) -> Callable:
"""Creates an expression for the quadrature of the right-hand-side of the
differential equation that is only dependent on the function values at the
nodes y_i. It takes the quadrature base points to evaluate the
base functions phi_i and their derivatives dphi_i to set up the expression.

Args:
rhs (Callable): Right hand side of the differential equation.
Needs to be a function of x, y, and dy
phi (List[np.poly1d]): List with all base polynomials phi. Mostly, base polynomials of Legendre type.
dphi (List[np.poly1d]): List with the corresponding derivatives of the base polynomials.
Returns:
Callable: Expression for the weighted quadratures of the right-hand-side for each node that is only
dependent on the unknown function values y_i. It can be called by a solving algorithm.
"""
# these will be defined during runtime
number_base_points = 2
base_points = np.array([0.00021132, 0.00078868])
weights = np.array([0.0005, 0.0005])
# evaluate base functions at quadrature base points
phi_eval = np.zeros(shape=(len(phi), number_base_points))
dphi_eval = np.zeros(shape=(len(dphi), number_base_points))
for index, phi_i in enumerate(phi):
phi_eval[index] = phi_i(base_points)
dphi_eval[index] = dphi[index](base_points)
# set up the expression for the quadrature only dependent on the unknown y_i values
def quad_expressions(y: np.array, model:GEKKO) -> np.array:
"""Variable expressions for the weighted quadratures of the integral of the weighted right-hand-side
for each node that are only dependend on y.
The integral is crucial for the Galerkin Method. int(phi_j * rhs(x, y, dy))dx.
It allows the gekko solver to embed a numpy array y into the quadrature
that is safed into the gekko model by using 'm.Array(m.Var)' for setting up the algebraic equation system.
Args:
y (np.array): Unknown function values that the gekko solver needs to find.
Returns:
np.array: Array with the weighted quadratures for each function and node as a function of y_i.
first index specifies the function, second index the node
"""
# TODO: Work with numpy arrays instead of lists to enhance speed.
y_tilde = []
dy_tilde = []
# set up the trail functions and their derivatives using the base functions which are already
# evaluated and the supposed y values at the nodes.
for eq_l in range(y.shape[0]):
y_tilde_eq = []
dy_tilde_eq = []
for bpoint_iq in range(number_base_points):
# iterates over quadrature base points
y_tilde_eq.append(np.dot(y[eq_l], phi_eval[:, bpoint_iq]))
dy_tilde_eq.append(np.dot(y[eq_l], dphi_eval[:, bpoint_iq]))
y_tilde.append(y_tilde_eq)
dy_tilde.append(dy_tilde_eq)
# evaluate the right-hand-sides with the base points and the supposed trial function values at the base points
rhsides_evaluated_bpoints = rhs(
base_points, np.array(y_tilde), np.array(dy_tilde)
)
# multiply the right-hand-side with the trial function of each node
integrals = np.zeros_like(y)
for eq_l, rhs_eq in enumerate(rhsides_evaluated_bpoints):
# iterate over the equations
for node_i in range(y.shape[1]):
# iterate over the weight function of each node
# calcualte the quadrature using the weights
rhs_w_func = rhs_eq * phi_eval[node_i]
dotp = model.sum(
[rhs_w_func[i] * weights[i] for i in range(number_base_points)]
)
# dotp = np.dot(
#     rhs_w_func, self.weights
# )
integrals[eq_l, node_i] = dotp  # using the dot product via list comprehension with gekko to allow the model to split long equations
return integrals
return quad_expressions

Setting up equations

phi = [
np.poly1d([ 2.e+06, -3.e+03, 1.e+00]),
np.poly1d([-4000000., 4000., -0.]),
np.poly1d([ 2.e+06, -1.e+03, 0.e+00])
]

dphi = [
np.poly1d([ 4.e+06, -3.e+03]),
np.poly1d([-8.e+06, 4.e+03]),
np.poly1d([ 4.e+06, -1.e+03])
]

diffusion_matrix = np.array(
[[ 2333.33333333, -2666.66666667, 333.33333333],
[-2666.66666667, 5333.33333333, -2666.66666667],
[ 333.33333333, -2666.66666667, 2333.33333333]]
)
y = m.Array(m.Var, (4, 3), value=x_bulk[0])

quad_expression = give_quad_expressions(rhs=rhs, phi=phi, dphi=dphi)
rhs_integrals = quad_expression(y=y, model=m)

for j in range(y.shape[0]):
# iterates over each differential equation
for i in range(y.shape[1]):
# iterates over each node
m.Equation(
np.dot(
y[j], (-1) * diffusion_matrix[:, i]
) == rhs_integrals[j, i]
)

for eq in m._equations:
print(eq)
m.solve()

Error
Exception: @error: Equation Definition
Equation without an equality (=) or inequality (>,<)
(-74900/((8.314)*(((((((v10)*(-0.12200771519999987))+((v11)*(0.6666554303999999)))+((v12)*(0.4553522848000001))))*(550)))))]))))*([((((((((v1)*(0.4553522848))+((v2)*(0.6666554304000001)))+((v3)*(-0.1220077152))))*(1500000)))^(0.6))
STOPPING...
</details>
# 答案1
**得分**: 2
我发现了问题：
出现这些换行符的原因是因为错误地使用了gekko内置函数，如`m.exp()`，`m.sin()`等。这些函数只能逐段评估，而不能像我试图做的那样一次性评估整个数组。
因此，从
```python
rhsides_evaluated_bpoints = rhs(
base_points, np.array(y_tilde), np.array(dy_tilde)
)

改为

rhs_eval = []
for i, bpoint in enumerate(base_points):
     rhs_eval.append(
          rhs(
            x=bpoint,
            y=y_tilde[:, i],
            dy=dy_tilde[:, i]
        )  
     )
rhs_eval = np.array(rhs_eval).transpose()

解决了问题。

rhsides_evaluated_bpoints = rhs(
base_points, np.array(y_tilde), np.array(dy_tilde)
)

rhs_eval = []
for i, bpoint in enumerate(base_points):
rhs_eval.append(
rhs(
x=bpoint,
y=y_tilde[:, i],
dy=dy_tilde[:, i]
)  
)
rhs_eval = np.array(rhs_eval).transpose()