# Python和MATLAB锂电铅蓄电化学微分模型和等效电路
1. 对比三种电化学颗粒模型：电化学的锂离子电池模型、单粒子模型和带电解质的单粒子模型。
2. 求解粒子域内边界通量与局部电流密度有关的扩散方程。
3. 扩展为两个相的负或正电极复合电极粒子模型。
4. 模拟四种耦合机制下活性物质损失情况。
5. 模拟锂离子电池三参数模型等效电路。

## :maple_leaf:电池解析模型

%%{init:{'gitGraph':{'mainBranchName':'电池'},'themeVariables':{'commitLabelBackground': 'None','commitLabelFontSize':'20px'}}}%%
  gitGraph
   commit id:"数学"
   commit id:"物理"
   commit id:"电化学"
   branch "汽车电池"
   commit id:"锂电池"
   branch "MATLAB电车(宝马i3)卸载电池再利用电气模型分析"
   branch "MATLAB和Python电车电池制造性能度量分析"
   commit id:"电化学指标"
   commit id:"热力学特性"
   commit id:"动力学特性"
   commit id:"活性材料"
   commit id:"定量分析"
   commit id:"锂化学计量"
   checkout "MATLAB电车(宝马i3)卸载电池再利用电气模型分析"
   commit id:"卸载电池"
   commit id:"电池容量"
   commit id:"电阻"
   commit id:"电气模型"
   checkout "汽车电池"
   commit id:"铅酸电池"
   branch "Python和MATLAB微分锂电铅蓄电化学模型和等效电路"
   commit id:"三种电化学颗粒模型"
   commit id:"扩散方程"
   commit id:"复合电极粒子模型"
   commit id:"活性物质损失"
   checkout "汽车电池"
   merge "MATLAB和Python电车电池制造性能度量分析"
   merge "MATLAB电车(宝马i3)卸载电池再利用电气模型分析"
   merge "Python和MATLAB微分锂电铅蓄电化学模型和等效电路"
   checkout "电池"
   commit id:"储能"
   branch "MATLAB(Octave)混电动力能耗评估"
   commit id:"储能元数据信息"
   commit id:"电流脉冲"
   commit id:"等效电路"
   commit id:"阻抗分析"
   checkout "电池"
   merge "汽车电池"
   merge "MATLAB(Octave)混电动力能耗评估"

## :cookie:语言内容分比

pie title 语言分比
 "Python":90
 "MATLAB":60
 "C/C++":30

pie title 内容分比
 "电化学、锂离子、铅酸、活性物质":60
 "微分模型":90
 "数学、常微分方程、偏微分方程、代数":80
 "物理学、电极、粒子、等效电路":50

## :grapes:Python和MATLAB符号微分
计算导数的常用方法有 3 种：

- 数值微分
- 符号微分
- 自动微分

数值微分依赖于导数的定义：
$$
f^{\prime}(x)=\lim _{h \rightarrow 0} \frac{f(x+h)-f(x)}{h}
$$
其中，你把一个非常小的 $h$ 和评估函数放在两个地方。这是最基本的公式，在实践中，人们使用其他公式，这些公式给出的估计误差更小。这种计算导数的方法主要适用于你不知道你的函数并且只能对其进行采样的情况。此外，对于高维函数，它需要大量的计算。

符号微分操纵数学表达式。如果你曾经使用过 matlab 或 mathematica，那么你会看到类似这样的内容

$$
\begin{aligned}
& \frac{d}{d x}\left(\frac{x^2 \cos (x-7)}{\sin (x)}\right)= \\
& x^2 \sin (7-x) \csc (x)+x^2(-\cos (7-x)) \cot (x) \csc (x)+2 x \cos (7-x) \csc (x)
\end{aligned}
$$
在这里，他们知道每个数学表达式的导数，并使用各种规则（乘积规则、链式规则）来计算最终的导数。然后他们简化表达式以获得最终表达式。

自动微分可操纵计算机程序块。微分器具有对程序的每个元素求导的规则。它还使用链式法则将复杂表达式分解为更简单的表达式。

### Python符号微分

导数规则要求一个项有 2 个子项，即二元运算（例如 *），以及一个一元运算（例如 sin）。

term.py:

from functools import reduce

class Term:

    def __init__(self, label: str or int, children=None):
        self.label = label
        self.children = children if children else []

    def __repr__(self):
        if self.children:
            return f"({self.label} " + reduce(lambda str1, str2: str1+" "+str2, (repr(child) for child in self.children)) + ")"
        else:
            return f"{self.label}"

    @classmethod
    def copy(cls, obj):
        """makes a copy of obj"""
        return cls(obj.label, [cls.copy(child) for child in obj.children])

symbols_table:

MUL = '*'
ADD = '+'
SUB = '-'
DIV = '/'
POW = '^'
LOG = 'log'
SIN = 'sin'
COS = 'cos'
TAN = 'tan'
OPEN_BR = '('
CLOSE_BR = ')'
VAR = 'x'

symbols = {LOG, SIN, COS, TAN, MUL, ADD, SUB, DIV, POW, OPEN_BR, CLOSE_BR, VAR}

from term import Term
from symbols_table import *
import numbers

def calculate_derivative(term1: Term, rules: dict) -> Term:
    if isinstance(term1.label, numbers.Number):
        return Term(0)
    elif term1.label == VAR:
        return Term(1)
    else:
        return rules[term1.label](term1, rules)


def mul_derivative(term1: Term, rules) -> Term:
    child1, child2 = term1.children  
    return Term(ADD, [Term(MUL, [calculate_derivative(child1, rules), child2]), Term(MUL, [child1, calculate_derivative(child2, rules)])])


def div_derivative(term1: Term, rules) -> Term:
    child1, child2 = term1.children  
    return Term(DIV, [Term(SUB, [Term(MUL, [calculate_derivative(child1, rules), child2]), Term(MUL, [child1, calculate_derivative(child2, rules)])]),
                      Term(POW, [child2, Term(2)])])


def sum_derivative(term1: Term, rules) -> Term:
    child1, child2 = term1.children  
    return Term(term1.label, [calculate_derivative(child1, rules), calculate_derivative(child2, rules)])


def power_derivative(term1: Term, rules) -> Term:
    child1, child2 = term1.children  
    return Term(ADD,
                [Term(MUL, [calculate_derivative(child2, rules), Term(LOG, [child1]), Term(POW, [child1, child2])]),
                 Term(MUL, [calculate_derivative(child1, rules), child2, Term(POW, [child1, Term(SUB, [child2, Term(1)])])])])


def sin_derivative(term1: Term, rules) -> Term:
    return Term(MUL, [Term(COS, [term1.children[0]]), calculate_derivative(term1.children[0], rules)])


def cos_derivative(term1: Term, rules) -> Term:
    return Term(MUL, [Term(-1), Term(SIN, [term1.children[0]]), calculate_derivative(term1.children[0], rules)])


def tan_derivative(term1: Term, rules) -> Term:
    arg = Term.copy(term1.children[0])
    return Term(ADD, [calculate_derivative(arg, rules), Term(MUL, [calculate_derivative(arg, rules), Term(POW, [term1, Term(2)])])])


def log_derivative(term1: Term, rules) -> Term:
    return Term(DIV, [calculate_derivative(term1.children[0], rules), term1.children[0]])


rules_for_differentiation = {MUL: mul_derivative, ADD: sum_derivative, SUB: sum_derivative,
                             COS: cos_derivative, SIN: sin_derivative, TAN: tan_derivative,
                             POW: power_derivative, DIV: div_derivative, LOG: log_derivative}


def derivative(term1):
    return calculate_derivative(term1, rules_for_differentiation)

### MATLAB符号微分

求函数 $y=f(x)^2 \frac{d f}{d x}$ 关于 $f(x)$ 的导数

syms f(x) y
y = f(x)^2*diff(f(x),x);
Dy = diff(y,f(x))

$\begin{aligned} & Dy = \\ & \qquad 2 f(x) \frac{\partial}{\partial x} f(x)\end{aligned}$

找到描述质量弹簧系统运动的欧拉-拉格朗日方程。定义系统的动能和势能。

syms x(t) m k
T = m/2*diff(x(t),t)^2;
V = k/2*x(t)^2;

定义拉格朗日量。

L = T - V

$\begin{aligned} & L = \\ & \frac{m\left(\frac{\partial}{\partial t} x(t)\right)^2}{2}-\frac{k x(t)^2}{2}\end{aligned}$

欧拉-拉格朗日方程如下

$0=\frac{d}{d t} \frac{\partial L(t, x, \dot{x})}{\partial \dot{x}}-\frac{\partial L(t, x, \dot{x})}{\partial x}$

计算项 $\partial L / \partial \dot{x}$

D1 = diff(L,diff(x(t),t))

$\begin{aligned} & D 1= \\ & m \frac{\partial}{\partial t} x(t)\end{aligned}$

计算第二项 $\partial L / \partial x$

D2 = diff(L,x)

$D 2( t )=-k x(t)$

找到质量弹簧系统的欧拉-拉格朗日运动方程。

diff(D1,t) - D2 == 0

$\begin{aligned}
& \operatorname{ans}(t)= \\
& m \frac{\partial^2}{\partial t^2} x(t)+k x(t)=0
\end{aligned}$

### :point_right: [亞圖跨際更多視角](https://viadean.notion.site/Python-MATLAB-1241ae7b9a32807e9b2afda4de53831a)