电路课程设计-模型构建

这个学期赶上了电路课的课改，作为计算机系的学生，要编程实现对任意结构线性电路的稳态、暂态分析和频率响应分析，硬核程度实在是让人觉得自己的程序设计课都白上了。

作为一个已经习惯了C++的退役选手，本着偷懒的想法没怎么考虑用Matlab和Python来做，结果偷懒不成反给自己找了不少麻烦，最终磕磕绊绊实现了一个还算能看的demo。考试周结束了，就随手总结一下自己这个差强人意的实践经历吧。

计划分成两篇来写，这一篇是关于数学模型和算法的，下一篇是关于程序设计的。不过我才疏学浅，如有错漏，敬请斧正。

以下是正文。

从图的角度对线性电路进行数学模型的建立和计算机辅助分析的设计，尝试解决线性电路在直流、交流激励下，电路的稳态分析问题和一类特殊线性动态电路的暂态分析问题。
思路是先实现稳态分析，再将连续的暂态分析离散化、近似线性化，视作稳态之间的转移过程。

线性电路的数学模型

名词约定

由于本文中同时设计电路和图论中一些相似的概念，区别如下:

• 所有的“结点”、“支路”、“回路”只代表电路中实际的结构;
• 所有的“点”、“边”、“环”只代表图论中抽象的元素。

线性电路的图描述

首先，对于课程涉及的线性电路而言，其导线的空间排布和实际长度都不对电路的任何物理量构成影响，关心的只有元件和元件之间的连接关系，这与图论中图的特点是一致的，所以选择用图来描述电路。

其次，虽然电路中的电流是具有方向的，但是此种方向只是人为规定，并不代表元件所建立的连接关系是单向的，所以选取的是电路的无向图表示而不是有向图表示。即，认为电路的图是无向图，所有的元件都对应无向边。

去掉具体元件的属性，线性电路中任何一个抽象的“元件”都可以用一条无向边表示，其必要的属性只有两个，即元件所连接的的两个结点，或图中边所连接的两个点。图论中对图所能进行的操作，例如遍历等，对电路也同样有效。

而对于实体的线性电路而言，每个元件还有一个参数，这个参数是决定该元件特性的物理量，所以将一个元件表示为

$$E=((m,n),c)$$

这是一个二元组，第一个元素$(m,n)$表示这个元件两端的结点，即其反映在图中的属性，第二个元素c表示元件的参数，即其实际的特性。

之所以不在元件的数学模型中加入端电压和电流，是因为端电压是完全可以由结点电压导出的，所以不作为元件的基本属性;而电流又可以通过 VCR 求出，不能用 VCR 求出的电流一定会被作为变量解出，所以电流也不是元件的基本属性。

稳态分析的求解

采用改进的结点电压法解决电路的稳态分析问题。（这部分可参见邱关源和罗先觉老师编写的《电路》教材）

已知电路的图的情况下，基本的结点电压方程表述为

$$Y\vec{u}=\vec{i}$$

其中电导(导纳)矩阵$Y$是实矩阵(表示电导)或复矩阵(表示导纳)，满足

$$Y(m,n)= \begin{cases} -\displaystyle{\sum_k^{Y_k=((m,n),c)}}Y_k, &m\neq n \\ -\displaystyle{\sum_k^{Y_k-((m,*),c)}}Y_k, &m=n \end{cases}$$

即第$m$行表示第$m$个结点的 KCL 方程，元素$Y(m,m)$是自导纳，其他元素都是互导纳。

变量$\vec{u}=(u_1,u_2,\cdots,u_k,\cdots)^{\mathrm{T}}$表示结点电压，列向量$\vec{i}=(i_1,i_2,\cdots,i_k,\cdots)^{\mathrm{T}}$表示各结点的注入电流。

对于结点电压法不允许无伴电流源的存在这个问题，作出如下的改进:

将流经无伴电压源的电流也作为变量加入方程组中，即将变量改为

$$\begin{bmatrix} \vec{u} \\ i_{VS} \end{bmatrix}$$

并且对于一个无伴电压源$VS=((m,n),U)$，增加一个方程

$$u_n-u_m=U$$

作为补充约束。

类似地，对于受控源，可以通过增加变量和增加方程的方式使之纳入到扩展的结点电压方程组中去。

最终形成的扩展结点电压方程组表示为

$$\begin{bmatrix} Y &S_1 \\ S_2 &S_3 \end{bmatrix} \begin{bmatrix} \vec{u} \\ \vec{p} \end{bmatrix} = \begin{bmatrix} \vec{i} \\ \vec{b} \end{bmatrix}$$

式中$S_1,S_2,S_3$都表示增加的方程带来的系数，而$\vec{p}$表示所有附加的变量，这些变量可能是流过电压源的电流，也可能是受控源的电压和电流。$\vec{b}$表示附加的方程右端的常数项。

如此，只需解该线性方程组，根据其解的情况可以给出线性电路的稳态分析结果：

• 如果方程组有解，那么电路是有稳态解的，并且可以求出任意结点的电压与任意支路的电流。
• 如果方程组无解，那么电路是没有稳态解的，电路中应当包含了使得某些元件不能工作的连接方式，例如短路的理想电压源或者开路的理想电流源。

暂态分析的求解

本文中尝试解决的暂态分析是对于线性动态元件而言的，即动态元件的特性（电容或电感）是常数，仅与元件自身的特性有关。

对于一个线性动态电路而言，如果满足如下的条件：

• 任何一个电容必定存在一个结点，该结点只连接有这一个电容
• 任何一个电感必定处于某一回路，该回路中只存在这一个电感

那么这个线性动态电路就可以根据如下的方法列出一阶的状态微分方程组：

1. 对于任一电容，找到其所连接的单电容结点，对这个结点列 KCL 方程，由于这个结点只连接了一个电容，个 KCL 方程必然是关于这个电容两端电压的一阶线性常微分方程
2. 对于任一电感，找到其所在的单电感回路，对这个回路列 KVL 方程，由于这个回路中只有一个电感，这个 KVL 方程必然是关于这个电感电流的一阶线性常微分方程

即，以电容电压和电感电流为状态变量，此类线性动态电路的状态方程一定可以表示为

$$\begin{bmatrix} \displaystyle{\frac{\mathrm{d}\vec{u}_C}{\mathrm{d}t}} \\\\ \displaystyle{\frac{\mathrm{d}\vec{i}_L}{\mathrm{d}t}} \end{bmatrix} =A \begin{bmatrix} \vec{u}_C \\ \vec{i}_L \end{bmatrix} +B\vec{y}$$

对于上述方程组，变量$\vec{y}$中可能含有激励源的参数或者电路中其他元件的电流、电压（中间变量），若将其整理为标准形式，即为教材中用激励作为输入所表示的状态方程。

但是，消去中间变量的过程是不易进行的。从数学推导的角度，这个过程需要对方程组进行反复加减调整，我没有能推导出形式化的方法；从程序实现的角度，激励源的电压（或电流）与元件的电压（或电流）在程序数据中并没有区别。

因此，我的思路是直接放弃消去中间变量，转而使用迭代的方式求微分方程组的数值解。

将时间离散化，取最小时间间隔为$\Delta t$，从而将上述的微分方程组写成差分方程组

$$\begin{bmatrix} \displaystyle\frac{\Delta\vec{u}_C}{\Delta t} \\\\ \displaystyle\frac{\Delta\vec{i}_L}{\Delta t} \end{bmatrix} =\frac{1}{\Delta t} \begin{bmatrix} \vec{u}_C(t_{n+1})-\vec{u}_C(t_n) \\ \vec{i}_L(t_{n+1})-\vec{i}_L(t_n) \end{bmatrix} =A \begin{bmatrix} \vec{u}_C(t_n) \\ \vec{i}_L(t_n) \end{bmatrix} +B\vec{y}(t_n)$$

根据极限理论，只要函数$\vec{u}_C(t)$和$\vec{i}_L(t)$是连续可导的，当$\Delta t\rarr 0$时，差商将收敛到微商，即可以用该差分方程组作为原微分方程组的一个近似，并且理论上可以使误差任意小。

令

$$\vec{x}(t)= \begin{bmatrix} \vec{u}_C(t) \\ \vec{i}_L(t) \end{bmatrix}$$

在时间离散化后，原先的时间轴变为时间序列$t_0=0,t_1,t_2,\cdots$，函数值也变为一个数列$\vec{x}(t_0),\vec{x}(t_1),\vec{x}(t_2),\cdots$，即状态的变化也是不连续的。那么，处于任一时刻$t_k$的电路就近似地成为稳态了。

由此，可以得到如下用于近似计算的性质：

在$t_k$时刻，电路近似稳态的依据是此时刻的状态量$\vec{x}(t_k)$近似不变。那么，若将$t_k$时刻的电路作为稳态电路近似处理，在用结点电压法求解这个近似的稳态电路时应满足

• 动态元件的状态量（即在原先不能突变的物理量）必须保持不变
• 原先能够突变的物理量则依然可以突变

根据元件的特性，$u_C$和$i_L$原先不能突变的，所以这两个物理量保持不变；$i_C$和$u_L$是原先可以突变的，所以这两个物理量依然可以突变。这样，在给定的时刻$t_k$，电容的特性完全类似于一个理想电压源，而电感的特性则类似于一个理想电流源了。

再观察方程组右端的形式，发现右端只含有$\vec{x}(t_n)$与$\vec{y}(t_n)$而不含有$\vec{x}(t_{n+1})$与$\vec{y}(t_{n+1})$。这就意味着只要知道了$\vec{x}(t_n)$与$\vec{y}(t_n)$的值，就可以直接递推计算出$\vec{x}(t_{n+1})$。而计算出$\vec{x}(t_{n+1})$后，根据电路在$t_{n+1}$时刻处于近似稳态，将电容视作电压源、电感视作电流源，又可以用结点电压法解出$\vec{y}(t_{n+1})$。

由此，只需确定初始条件$\vec{x}(t_0)$，即可计算出电路在之后任意时刻的状态变量$\vec{x}(t_k)$，将计算出的（离散）状态变量序列作分段线性插值所得的连续曲线，就是这个微分方程组的一组数值解曲线。

算法分析与设计

稳态分析的求解

稳态分析的求解就是填写结点电压方程系数矩阵和求解系数矩阵的过程。

填写矩阵的方式与教材基本相同，另需提前统计附加变量的个数，并根据附加的方程填入合适的位置。

线性方程组的求解则用高斯-约旦消元法，该方法与高斯消元法的不同在于：此方法会将系数矩阵消元至对角矩阵而不是上三角矩阵，相对来说回代次数较少，精度更高一些。

单电容节点的判定

直接统计所有的节点上链接的电容个数，就可以判定出单电容节点。

单电感回路的判定

单电感回路分为两种情况：

1. 回路含有 2 个结点的，对应于图上的二元环
2. 回路含有多于 2 个结点的，对应于图上的多元环

对于第一种情况，只需要枚举所有点对，判定它们之间是否存在一个环，对于第二种情况，则用深度优先遍历(DFS)的方式找到多元环：

1. 从一个点开始遍历，每找到一条边
2. 如果所邻接的点没有被访问过，就将邻接的点入栈，标记这个点已经被访问，并访问这个点
3. 如果所邻接的点已经被访问过，并且该点在栈中，又不是当前点在搜索树中的父节点，则说明已找到一个多元环，依次退栈直至栈顶的点就是当前的邻接点为止，这些出栈的点就构成了一个环
4. 检查所有的边后，该点退栈，回溯

这个在无向图中寻找多元简单环的方法基于一个结论：无向图的深度优先搜索树中一定只含有树边和后向边，不含有前向边和横叉边，于是简单环一定是由一些树边和一条后向边组成的。

暂态分析的计算

并不显式地填写差分方程的系数矩阵$A$和$B$，而是对每个动态元件都维护一个表，包含所有与之相关的元件（及这些元件的电压、电流信息），每次从$\vec{x}(t_n)$计算$\vec{x}(t_{n+1})$时就提取这些原件的信息，算出$\vec{x}(t_{n+1})$后再更新所有元件的电压和电流信息。

总结和讨论

收获与感想

我对软件的兴趣要大于对硬件的兴趣，所以电路实在不算是我比较喜欢的一类课程，只能说是并不讨厌。（当然，也许是我所学尚浅，不能领会其中深意，如有指教，欢迎交流）

对于这学期的任务，我抱有一点小小的奢望：能去“设计一个用数值计算的方法解决带有工程（物理）背景问题的软件”、培养自己“拿到一些工程知识之后，为工程问题设计合适的软件和数据结构”的能力。

我进入计算机系之前的身份是OI选手，考虑问题一般只需要关心计算复杂度（甚至往往只关心时间复杂度就行了），缺少软件整体设计的训练和思考。这次完成任务是我第一次做出一些不同的尝试：让自己的程序更模块化，有更好的可重用、可扩展性。

在第一次任务的成果中，我将解线性方程组的部分完全独立出来，并且初步设计了数据、元件、电路这三个抽象层次（这些内容会在下一篇中谈到）。后来的第二次和第三次任务基本上直接复用了求解方程组的部分，也沿用了元件的封装结构、接口，又进一步将电路这个层次比较完整地抽象出来。并且，对简单的线性电阻电路增加暂态分析功能和相量分析功能的时候，也在基本不改动原先体系、尽量利用原有功能的前提下完成了扩展工作。在整个任务过程中，我尝试着将所学的OO思想应用到编程实践中，也体会到了多态、封装这两个OO特性带来的独特优势，有意识地将软件视作一个系统去协调各个层次的关系，这些都引发过我不少思考，让我获益良多。

另外，十分感谢负责专题实验的老师对我的指导，虽然我们在具体执行上没有过交流，但是老师以学长的身份对我的学习生活给出了诚恳的建议。处于保护隐私的需要，老师的身份和交流的经历不再详谈。

不足之处

这次实验设计比较有难度的是动态电路的暂态分析问题，而我求解暂态分析的思路是基于教材 7-10 一节所介绍的状态方程的。

最初我希望能够列出状态方程组，再用比较经典的方法（比如龙格-库塔法等等）去求数值解，这样精度比较好。但是在消去中间变量的过程中遇到了困难，又因为时间比较紧张，没有能更多地参考和学习（电路）网络图论的知识，以决定转换思路，用一阶线性近似，避开消去中间变量的问题。

我应用状态方程的思路决定了设计的程序无法求解不能列出一阶状态微分方程组的电路（比如，哪怕是最简单的 2 个并联电容）。如果要补上这个不足，就需要彻底更换方法。这不得不说是一个遗憾。