学习数值分析与电路分析后的练手程序，较粗糙。 没有 加入更多电器元件的数学模型，如二极管、三极管、放大器、数字器件等等。
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 实现。

输入与输出

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

电路结构

信号源

仿真时间序列

初始值

输出

数据处理

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

初始化

• 根据输入变量的个数判断所描述电路的类型
  • 根据电路求解类型，恰当设置其余变量的值
  • 没有用到的变量赋予空，如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. 加入更多电器元件的数学模型，如二极管、三极管、放大器、数字器件等等

最后

以上就是舒适香烟为你收集整理的【matlab】数值计算实现电路仿真搞清几个概念思路具体实行编程实现拓展的全部内容，希望文章能够帮你解决【matlab】数值计算实现电路仿真搞清几个概念思路具体实行编程实现拓展所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。