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

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

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()