【matlab】数值计算实现电路仿真

学习数值分析与电路分析后的练手程序，较粗糙。 没有 加入更多电器元件的数学模型，如二极管、三极管、放大器、数字器件等等。
m脚本在：gittee仓库；

---

编程实现电路仿真，使用Matlab 。
背景：
电网络分析课程中，学习了电路的矩阵形式，和数值求解方法；数值分析课程了解了数值计算的概念和一些数值算法思想。于是我想，如何编程实现电路的仿真呢？

---

搞清几个概念

总是理不清思路：编程计算电路？理论分析电路？怎么转换呢？有什么不一样吗？得到各参数变量与待求量之间的关系？得到不同电路参数下的解？？？
所以首先要明确模糊的概念，以便理清思路。

数值计算

利用计算机求解数学问题，称为数值计算。计算机可以大量、快速的运算，但是它是基于硬件结构的机械式计算，不同于人的计算（更高级、更灵活、更包容、更连续）。

需要考虑到计算机计算与编程的特点，设计恰当的数学模型、合适的数值计算方法和形式。

求解

求解即已知题目求解答案。有得到精确解析解的数值计算方法，也有迭代近似解的数值方法。计算只是个过程，重要的是，题目是什么，想要得到什么样的结果？

需要确定什么是已知的，什么是未知的。再理清求解思路：先分析，列出已知条件，然后计算出结果，再进一步得到更多结果，这才是完整的逻辑。这个逻辑是理论（课程）中学到的，这说明理论与求解是相辅相成的，理论上的分析虽然不需要数值计算，但没有算例验证理论，可信度也不高；同时，算例数据可以直观的反映出问题的特点，给分析带来参考，有助于得出结论。

仿真

计算机仿真：一般以时间为变量，将仿真过程分成极短的片段，求解出每一时刻的状态近似作为对应片段内的状态。从而近似得到整个过程中任意时刻的状态。

而求解得到的是一个问题的答案（一个电路状态信息），而仿真需要每个时间片段的答案（每个时刻电路的状态信息），这意味着“连续求解”一个“连续的题目”，需要将前一个问题的答案作为已知去求解当前问题。

一次仿真，得到的往往是大量的求解，仿真的结果应当是连续的（连续的时间序列）、大量的解的集合（按照时间序列依次求解的解序列）。

---

思路

清楚上面的概念后，我理清了思路：

1. 理论分析电路，编写能求解“某时刻电路”的程序（单次求解）；
2. 加入相邻时刻间的变量更新，在时间序列的每个时刻求解，实现电路仿真（连续多次求解）；
3. 改变电路参数，得到相应结果，解析法验证（计算算例，通过理论计算验证）。

---

具体实行

这是最难的一部分：概念上，混淆了求解与仿真；思想上，纠结于解析法，而忽略数值方法；实现上，多方面知识结合。

描述瞬时电路

利用图论的概念描述电网络，形成网络线图（G图）：

• 抽象出节点和支路
• 每一元件为一条支路
• 各支路标有方向箭头（自指定方向）

利用电网络理论的知识，可得到电路的关联矩阵 $\mathbf{A}$，用于描述电路的拓扑结构。

• A = { a i j } {\boldsymbol A} = \{a_{ij}\} A={aij​}
• $a_{ij}=1$：支路 j 箭头指向节点 i
• $a_{ij}=-1$：支路 j 箭头背离节点 i
• $a_{ij}=0$：支路 j 不与节点 i 关联

若有N+1个节点，M条支路，a个电导（电阻），b个电压源（包括电容），c个电流源（包括电感），
并考虑到极短时间内各参数不变，有以下矩阵：

• ${\mathbf{U}}_{\mathrm{s}}$：支路源电压向量（已知）
  • ${\mathbf{U}}_{\mathrm{s}}=[u_{\mathrm{s}1},u_{\mathrm{s}2},\cdots ,u_{\mathrm{s}\mathrm{M}}{]}^{\mathrm{T}}$
  • $u_{\mathrm{s}k}$ ：第 k 条支路的电压源电压，符号以G图对应支路的参考方向为正；k = 1 ~ M
  • 电容视为电压源，其值为初始电压
  • 无电压源或电容则为0
• ${\mathbf{I}}_{\mathrm{s}}$：支路源电流向量（已知）
  • ${\mathbf{I}}_{\mathrm{s}}=[i_{\mathrm{s}1},i_{\mathrm{s}2},\cdots ,i_{\mathrm{s}\mathrm{M}}{]}^{\mathrm{T}}$
  • $i_{\mathrm{s}k}$：第 k 条支路的电流源电流，符号以G图对应支路的参考方向为正；k = 1 ~ M
  • 电感视为电流源，其值为初始电流
  • 无电流源或电感则为0
• $\mathbf{Y}$：支路电导矩阵（已知）
  • $\mathbf{Y}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}[y_1{y}_2\cdots y_{\mathrm{M}}]$
  • $y_k$：第 k 条支路的电导；k = 1 ~ M
  • 电阻转化为电导，无电导则为0
  • 好像含受控源的电路的受控关系可以用电导矩阵表示出来（没研究）
• $\mathbf{U}$：支路电压向量（未知）
  • $\mathbf{U}=[u_1,u_2,\cdots ,u_{\mathrm{M}}{]}^{\mathrm{T}}$
  • $u_k$：第 k 条支路的电压，符号以G图对应支路的参考方向为正；k = 1 ~ M
• $\mathbf{I}$：支路电流向量（未知）
  • $\mathbf{I}=[i_1,i_2,\cdots ,i_{\mathrm{M}}{]}^{\mathrm{T}}$
  • $i_k$：第 k 条支路的电流，符号以G图对应支路的参考方向为正；k = 1 ~ M
• ${\mathbf{U}}_{\mathrm{n}}$：节点电压向量（未知）
  • ${\mathbf{U}}_{\mathrm{n}}=[u_{\mathrm{n}1},u_{\mathrm{n}2},\cdots ,u_{\mathrm{n}\mathrm{N}}{]}^{\mathrm{T}}$
  • $u_{\mathrm{n}k}$：第 k 个节点的电压；k = 0 ~ N
  • $u_{\mathrm{n}0}$：0号节点为参考节点

求解瞬时电路

将电容器和电感器分别看作电压源和电流源，“瞬间”（极短的时间，不是时刻）电路近似为一个只有电源和电阻组成的直流电路。

我只想到了网孔电流法和节点电压法，相比之下，节点电压法不需要“确定网孔”环节，故使用节点电压法。

由电路的 KCL 可知：
$$\mathbf{A}\mathbf{I}=0\text{\hspace{1em}(1)}$$
支路电压可用节点电压表示：
$$\mathbf{U}={\mathbf{A}}^{\mathrm{T}}{\mathbf{U}}_{\mathrm{n}}\text{\hspace{1em}(2)}$$
支路电流可以分为电导支路电流、电压源支路电流、电流源支路电流三个部分：
$$\mathbf{I}=\mathbf{Y}\mathbf{U}+{\mathbf{I}}_{\mathrm{s}}+{\mathbf{I}}_{\mathrm{u}}\text{\hspace{1em}(3)}$$
其中${\mathbf{I}}_{\mathrm{u}}$是除电压源支路以外的支路电流均为0的支路电流向量。

联立(1)(2)(3)，消去的矩阵$\mathbf{U}$和$\mathbf{I}$，得
$${\mathbf{A}\mathbf{Y}\mathbf{A}}^{\mathrm{T}}{\mathbf{U}}_{\mathrm{n}}+{\mathbf{A}\mathbf{I}}_{\mathrm{u}}+{\mathbf{A}\mathbf{I}}_{\mathrm{s}}=0\text{\hspace{1em}(4)}$$

在(4)表示的方程组中，${\mathbf{U}}_{\mathrm{n}}$与${\mathbf{I}}_{\mathrm{u}}$是待求量，有 N + b 个；独立方程数是 N 个（$\mathbf{A}$行数）。
显然，${\mathbf{U}}_{\mathrm{s}}$这个已知量还没有用到。

支路电压也可分为电导支路电压、电压源支路电压、电流源安支路电压三个部分：
$$\mathbf{U}={\mathbf{U}}_{\mathrm{y}}+{\mathbf{U}}_{\mathrm{s}}+{\mathbf{U}}_{\mathrm{i}}\text{\hspace{1em}(5)}$$
其中${\mathbf{U}}_{\mathrm{y}}$和${\mathbf{U}}_{\mathrm{i}}$分别是除电导支路以外的支路电压均为0的支路电压向量，和除电流源支路以外的支路电压均为0的之路电压向量。

由(2)和(5)可得，
$${\mathbf{A}}^{\mathrm{T}}{\mathbf{U}}_{\mathrm{n}}={\mathbf{U}}_{\mathrm{y}}+{\mathbf{U}}_{\mathrm{s}}+{\mathbf{U}}_{\mathrm{i}}$$
取出该 方程组 中，等号右侧是已知的$u_{\mathrm{s}}$的方程，表示为：
$$\mathbf{B}{\mathbf{U}}_{\mathrm{n}}={\mathbf{U}}_{\mathrm{s}}^{\prime }\text{\hspace{1em}(6)}$$
其中，$\mathbf{B}$和${\mathbf{U}}_{\mathrm{s}}^{\prime }$分别是${\mathbf{A}}^{\mathrm{T}}$和${\mathbf{U}}_{\mathrm{s}}$删除${\mathbf{U}}_{\mathrm{s}}$中0所在行得到的。
(6)所表示的方程组中，未知量是${\mathbf{U}}_{\mathrm{n}}$，独立方程数为 b ，即电压源支路的个数。

将(4)和(6)组成方程组，
$$\{\begin{array}{c}{\mathbf{A}\mathbf{Y}\mathbf{A}}^{\mathrm{T}}{\mathbf{U}}_{\mathrm{n}}+{\mathbf{A}\mathbf{I}}_{\mathrm{u}}=-{\mathbf{A}\mathbf{I}}_{\mathrm{s}}\\ \mathbf{B}{\mathbf{U}}_{\mathrm{n}}={\mathbf{U}}_{\mathrm{s}}^{\prime }\end{array}$$
其矩阵形式为：
$$\left[\begin{array}{cc}{\mathbf{A}\mathbf{Y}\mathbf{A}}^{\mathrm{T}}& \mathbf{A}\\ \mathbf{B}& 0\end{array}\right]\left[\begin{array}{c}{\mathbf{U}}_{\mathrm{n}}\\ {\mathbf{I}}_{\mathrm{u}}\end{array}\right]=\left[\begin{array}{c}-{AI}_{\mathrm{s}}\\ {\mathbf{U}}_{\mathrm{s}}^{\prime }\end{array}\right]\text{\hspace{1em}(7)}$$
当当当！ 该方程组有 N + b 个未知量，N + b 个独立方程，可解出${\mathbf{U}}_{\mathrm{n}}$和${\mathbf{I}}_{\mathrm{u}}$。

还可以进一步精简${\mathbf{A}\mathbf{I}}_{\mathrm{u}}$部分的形式：
根据${\mathbf{I}}_{\mathrm{u}}$中0行的位置，删去${\mathbf{A}\mathbf{I}}_{\mathrm{u}}$部分中系数矩阵$\mathbf{A}$无用的列，和${\mathbf{I}}_{\mathrm{u}}$中无用的行，
表示为${\mathbf{D}\mathbf{I}}_{\mathrm{u}}^{\prime }$，则(7)可写为
$$\left[\begin{array}{cc}{\mathbf{A}\mathbf{Y}\mathbf{A}}^{\mathrm{T}}& \mathbf{D}\\ \mathbf{B}& 0\end{array}\right]\left[\begin{array}{c}{\mathbf{U}}_{\mathrm{n}}\\ {\mathbf{I}}_{\mathrm{u}}^{\prime }\end{array}\right]=\left[\begin{array}{c}-{AI}_{\mathrm{s}}\\ {\mathbf{U}}_{\mathrm{s}}^{\prime }\end{array}\right]\text{\hspace{1em}(8)}$$

将它们分别回代到(2)和(3)则得到电路的此时的所有信息。

更新动态元件初始值

对于含动态元件（一般为电容、电感）的电路，每个瞬间里动态元件的电压和电流都在改变。
所以在每个瞬间电路求解完后，要利用求得的解对动态元件的初始值更新，作为下一瞬间动态元件的初值状态。

易知电容和电感的微分方程可写为：
$$\begin{gathered} \mathrm{d}{u}_{\mathrm{C}}=\frac{1}{C}{i}_{\mathrm{C}}\mathrm{d}t \\ \mathrm{d}{i}_{\mathrm{L}}=\frac{1}{L}{u}_{\mathrm{L}}\mathrm{d}t \end{gathered}$$
其中$\mathrm{d}t$可近似看作是一个“瞬间”的时间，在计算中称为步长，设为$\Delta t$，即$\mathrm{d}t=\Delta t$。
则得到电容电压变化量$\Delta u_{\mathrm{C}}$和电感电流变化量$\Delta i_{\mathrm{L}}$的表达式：
$$\Delta u_{\mathrm{C}}=\frac{1}{C}{i}_{\mathrm{C}}\Delta t\text{\hspace{1em}(8)}$$

$$\Delta i_{\mathrm{L}}=\frac{1}{L}{u}_{\mathrm{L}}\Delta t\text{\hspace{1em}(9)}$$

其中$i_{\mathrm{C}}$和$u_{\mathrm{L}}$分别是电容支路的电流和电感支路的电压，可从$\mathbf{I}$与$\mathbf{U}$中得到。

每次求解后，对每个电容的初始电压根据(8)式更新，每个电感的初始电流根据(9)式更新：
$$u_{\mathrm{C}}=u_{\mathrm{C}}+\Delta u_{\mathrm{C}}\text{\hspace{1em}(10)}$$

$$i_{\mathrm{L}}=i_{\mathrm{L}}+\Delta i_{\mathrm{L}}\text{\hspace{1em}(11)}$$

---

编程实现

使用 matlab 实现。

输入与输出

仿真电路需要输入：电路结构模型、信号源（电源）信号、仿真时间、动态元件初始值；
输出：电路的节点电压、支路电压、支路电流矩阵，包含历史记录。

电路结构

变量	规模	说明	意义
$\mathbf{A}$	$n\times b$	节点支路的关联情况 行数 = 节点数，列数 = 支路数	电路线图的关联矩阵
$\mathbf{Y}$	$b$	支路电导分布在对角线上 不可为空向量	支路电导向量
$\mathbf{C}$	$b$	支路上的电容值，单位 F 非电容支路值为 0	支路电容向量
$\mathbf{L}$	$b$	支路上的电感值，单位 H 非电感支路值为 0	支路电感向量

信号源

变量	规模	说明	意义
${\mathbf{U}}_{\mathrm{s}\mathrm{s}}$	$b$	支路恒压源电压值，单位 V 非恒压源支路值为 0 不可为空向量	支路恒压源电压向量
${\mathbf{I}}_{\mathrm{s}\mathrm{s}}$	$b$	支路恒流源电流值，单位 V 非恒流源支路值为 0 不可为空向量	支路恒流源电压向量
${\mathbf{U}}_{\mathrm{d}\mathrm{s}}$	$uds\times (1+T)$	第一列是索引 其余列是电压源随时间的输入数据 单位 V $uds$：变化的电压源数 $T$：时间步数	指定变化电压源的 位置和时变信号
${\mathbf{I}}_{\mathrm{d}\mathrm{s}}$	$ids\times (1+T)$	第一列是索引 其余列是电流源随时间的输入数据 单位 A $uds$：变化的电流源数 $T$：时间步数	指定变化电流源的 位置和时变信号

仿真时间序列

变量	规模	说明	意义
${\mathbf{T}}_{\mathrm{s}\mathrm{e}\mathrm{q}}$	$T$	时刻组成的序列，单位 s	仿真的时间序列

初始值

变量	规模	说明	意义
${\mathbf{U}}_c$	$c$	仿真开始前，各电容的初始值	电容初始值的向量
${\mathbf{I}}_l$	$l$	仿真开始前，各电感的初始值	电感初始值的向量

输出

变量	规模	说明	意义
${\mathbf{U}}_{\mathrm{n}}$	$n\times T$	每一时刻的节点电压向量组成的矩阵	节点电压矩阵
$\mathbf{U}$	$b\times T$	每一时刻的支路电压向量组成的矩阵	支路电压矩阵
$\mathbf{I}$	$b\times T$	每一时刻的支路电流向量组成的矩阵	支路电流矩阵

数据处理

处理输入信息，生成必要中间变量，构建求解矩阵。

初始化

• 根据输入变量的个数判断所描述电路的类型
  • 根据电路求解类型，恰当设置其余变量的值
  • 没有用到的变量赋予空，如U_ds = []
• 输入变量形式转化
  • 转化为向量形式，如U_ss = U_ss(:)等
  • 转化为矩阵形式，Y = diag(Y)
• 提取出${\mathbf{T}}_{\mathrm{s}\mathrm{e}\mathrm{q}}$的长度与步长
• 确定$\mathbf{U}$和$\mathbf{I}$矩阵的规模，即是输出矩阵，又作为相关变量的存储矩阵
  • 将已知的支路电压和支路电流填入其中
  • 在填入前，应判断变量是否存在，或填入空的变量是否影响计算，再决定是否填入
  • ${\mathbf{U}}_{\mathrm{s}\mathrm{s}}$是支路恒压源电压向量，长度 = 支路数，不随时间变化
    • 可用来创建$\mathbf{U}$ ：U = repmat(U_ss, 1, T)
    • 同时，要创建恒压源支路的索引：U_ss_idx = find(U_ss)
  • ${\mathbf{U}}_{\mathrm{d}\mathrm{s}}$是变压源向量，第一列是支路位置索引，其余是对应时刻数值
    • 提取变压源索引：U_ds_idx = U_ds(:, 1)
    • 放入$\mathbf{U}$：U(U_ds_idx,:) = U_ds(:, 2:end)
  • ${\mathbf{U}}_{\mathrm{c}}$是电容初始值向量
    • 可通过支路电容向量$\mathbf{C}$找到其支路位置索引：C_idx = find(C)
    • 放入$\mathbf{U}$：U(C_idx, 1) = U_c
  • ${\mathbf{I}}_{\mathrm{s}\mathrm{s}}$是支路恒流源电流向量，长度 = 支路数，不随时间变化
    • 可用来创建$\mathbf{I}$ ：I = repmat(I_ss, 1, T)
    • 同时，要创建恒压源支路的索引：I_ss_idx = find(I_ss)
  • ${\mathbf{I}}_{\mathrm{d}\mathrm{s}}$是变压源向量，第一列是支路位置索引，其余是对应时刻数值
    • 提取变压源索引：I_ds_idx = I_ds(:, 1)
    • 放入$\mathbf{I}$：I(I_ds_idx,:) = I_ds(:, 2:end)
  • ${\mathbf{I}}_{\mathrm{l}}$是电感初始值向量
    • 可通过支路电感向量$\mathbf{L}$找到其支路位置索引：L_idx = find(L)
    • 放入$\mathbf{I}$：I(L_idx, 1) = I_l
• 除了上述索引，还有构建求解方程组要用到的索引
  • 电压源和电容支路的索引：U_sc_idx = [U_ss_idx; C_idx; U_ds_idx]
  • 电流源和电感支路的索引：I_sl_idx = [I_ss_idx; L_idx; I_ds_idx]

构建方程组的系数矩阵

系数矩阵固定不变，方程组(8)的系数矩阵为：
$$\left[\begin{array}{cc}{\mathbf{A}\mathbf{Y}\mathbf{A}}^{\mathrm{T}}& \mathbf{D}\\ \mathbf{B}& 0\end{array}\right]$$
实现方式为：

AT = A'; % A 的转置矩阵
Y = diag(Y); % 电导向量转化为对角矩阵
B = AT(U_sc_idx,:); % 已知的节点电压关系矩阵
D = A(:,U_sc_idx); % KCL 中电压源和电容支路电流的系数矩阵
R = [A*Y*AT D; B zeros(length(B(:,1)),length(D(1,:)))]; % 构建求解方程组的矩阵

迭代求解与记录

使用for循环遍历时间序列。每一时刻（每一步）都进行一次求解，将结果填入${\mathbf{U}}_{\mathrm{n}}\mathbf{U}\mathbf{I}$矩阵的相应列。

方程的常数项

方程组的常数项向量为：
$$\left[\begin{array}{c}-{AI}_{\mathrm{s}}\\ {\mathbf{U}}_{\mathrm{s}}^{\prime }\end{array}\right]$$
常数项向量中的变化源、电容电压、电感电流是随时间变化的，
每次求解前都需要读取当前时刻(t时刻)的值：

    I_sl = zeros(b, 1);           % I_sl 电流源和电感支路的电流向量
    I_sl(I_sl_idx) = I(I_sl_idx,t); % 更新常数项
    U_sc = U(U_sc_idx,t);           %

I_sl即是${\mathbf{I}}_{\mathrm{s}}$，长度为 $b$，除电流源和电感支路的电流以外，其他支路均为0。
U_sc即是${\mathbf{U}}_{\mathrm{s}}^{\prime }$，长度为电压源和电容个数之和，可以为空；电压源和电容的电压向量。

求解和记录

求解方程组得到的向量是：
$$\left[\begin{array}{c}{\mathbf{U}}_{\mathrm{n}}\\ {\mathbf{I}}_{\mathrm{u}}^{\prime }\end{array}\right]$$
${\mathbf{U}}_{\mathrm{n}}$是此刻的节点电压向量，${\mathbf{I}}_{\mathrm{u}}^{\prime }$是此刻电压源和电容的电流；
根据方程组(2)得到此刻的$\mathbf{U}$，再根据(3)得到此刻的$\mathbf{I}$；
最后通过公式(10)、(11)分别更新电容的电压值、电感的电流值，作为下一时刻的初始值。

实现：

    temp = R\[-A*I_sl; U_sc];       % 求解方程
        
    U_n(:,t) = temp(1:sizeA(1));    % 得到 t 时刻的 U_n
    U(:,t) = AT*U_n(:,t);           % 得到 t 时刻的 U

    temp(1:sizeA(1)) = [];      % 删除结果中 U_n 的部分
    I_sc = temp;                % 得到 I_sc
    
    I(U_sc_idx,t) = I_sc;        % 将 t 时刻的 I_sc 添加到 I
    I(:,t) = I(:,t) + Y*U(:,t);  % 得到 t 时刻的 I
        
    if t < length(T_seq)	% 最后一次不需要计算电容或电感的值
        a = T_seq(t+1) - T_seq(t);			   % 计算步长
        U_c = U_c + a.*I(C_idx,t)./C(C_idx);    % 更新 U_c   I(C_idx,t)即是电容的电流向量
        I_l = I_l + a.*U(L_idx,t)./L(L_idx);    % 更新 I_l   U(L_idx,t)即是电感的电压向量
        U(C_idx, t+1) = U_c;
        I(L_idx, t+1) = I_l;
    end

其中，I_sc即是${\mathbf{I}}_{\mathrm{u}}^{\prime }$，电压源和电容的电流向量；
完成最后一次求解后，不再需要计算下一时刻电容或电感的值，使用if实现。

---

拓展

1. 含受控源的电路应该可以用此方法计算
2. 对于电路拓扑结构统一的电路（如，分布参数电路），可以快速灵活改变 电路结构 信息，分析拓扑结构相同而节点、支路数量不同对电路的影响
3. 加入更多电器元件的数学模型，如二极管、三极管、放大器、数字器件等等