x_Iij_square=sdpvar(32,N,'full');%电流平方 x_ui_square=sdpvar(33,N,'full');%电压平方 x_pij=sdpvar(32,N,'full');%线路有功功率 x_qij=sdpvar(32,N,'full');%线路无功功率 x_Pdg=sdpvar(3,N,'full');%dg功率 u_Pdg=binvar(3,N,'full');%dg出力状态 x_Pil=sdpvar(2,N,'full');%切负荷功率 u_Pil=binvar(2,N,'full');%切负荷状态 x0_ug=binvar(3,N,'full');%机组启动 x0_vg=binvar(3,N,'full');%机组停止 x_R=sdpvar(1,N,'full');%备用负荷 %% 定义约束条件 Constraints=[]; P_dg=sdpvar(33,24,'full');%将dg出力情况放入对应的节点 for i=1:33 if i==7 Constraints=[Constraints,P_dg(i,:)==x_Pdg(1,:)]; elseif i==24 Constraints=[Constraints,P_dg(i,:)==x_Pdg(2,:)]; elseif i==25 Constraints=[Constraints,P_dg(i,:)==x_Pdg(3,:)]; else Constraints=[Constraints,P_dg(i,:)==0]; end end P_il=sdpvar(33,24,'full');%将IL合同约束放入对应的节点 for i=1:33 if i==8 Constraints=[Constraints,P_il(i,:)==x_Pil(1,:)]; elseif i==25 Constraints=[Constraints,P_il(i,:)==x_Pil(2,:)]; else Constraints=[Constraints,P_il(i,:)==0]; end end Constraints=[Constraints,x_ui_square(1,:)==12.66^2]; %平衡节点每小时电压平方 Constraints=[Constraints,x_Iij_square>=0]; % Constraints=[Constraints,x_pij(1,:)>=0]; %% 支路欧姆定律 for r=1:32 Constraints=[Constraints,x_ui_square(Branch(r,2),:)-... x_ui_square(Branch(r,3),:)+(r_ij(r)^2+x_ij(r)^2)*x_Iij_square(r,:)-... 2*(r_ij(r)*x_pij(r,:)+x_ij(r)*x_qij(r,:))==0]; end %dg爬坡约束 for t=1:23 for i=1:3 Constraints=[Constraints,x_Pdg(i,t+1)-x_Pdg(i,t)<=vup]; Constraints=[Constraints,x_Pdg(i,t)-x_Pdg(i,t+1)<=vdn]; end end %IL合同约束 Constraints=[Constraints,0<=x_Pil<=Pilmax.*u_Pil]; %dg最小启停时间约束 for i=1:3 Constraints=[Constraints,consequtiveON([onoffhistory(i,:) u_Pdg(i,:)],Tmup)]; Constraints=[Constraints,consequtiveON(1-[onoffhistory(i,:) u_Pdg(i,:)],Tmd)]; end for t=2:N%公式7这段代码什么意思

<think>好的，我现在要解决的问题是，如何根据用户提供的MATLAB插值代码，将变量x_Iij_square进行类似的插值处理。首先，我需要仔细理解用户提供的代码片段，然后分析他们的变量结构和插值方法，最后应用到新的变量上。 首先，用户给出的代码中，x被定义为1到N的整数序列，间隔为1。然后x1被定义为从1开始，以1/4为间隔，直到N的序列。这意味着x1的长度是原来的四倍，比如原来的N个点，现在每个间隔被分成四个小段，所以总共有4*(N-1)+1个点？或者需要确认具体长度。比如，如果原x是1:1:N，即步长1，那么x1是1:0.25:N，例如当N=5时，x1会有17个点（1, 1.25, 1.5,...,5）。这相当于将每个原始间隔分成了四份，时间分辨率提高了四倍。 接下来是循环部分，nk从1到102，循环处理p_Solar1和p_Load1。这里，p_Solar_1和p_Load_1似乎是原始数据，维度可能是102行，每行有N个数据点。然后通过interp1函数，使用'spline'样条插值方法，将每个nk对应的行数据从x插值到x1上，结果存储在p_Solar1(nk,:)和p_Load1(nk,:)中，所以每个nk对应的行会有更多的列，对应于x1的长度。 现在的问题是，用户有一个变量x_Iij_square，它是通过sdpvar定义的，维度是32行N列的矩阵。用户希望按照同样的插值方法对这个变量进行插值处理。这里可能需要明确的是，x_Iij_square是一个优化变量还是已经存在的数据？根据sdpvar的用法，这应该是YALMIP中的决策变量，用于优化问题建模。但插值通常是对已知数据进行处理，所以这里可能需要确认是否在优化过程中需要生成插值后的变量，或者在优化问题中使用插值后的变量？ 假设用户想要对x_Iij_square的每一行进行插值，类似于p_Solar1的处理方式。即，原数据是32行，每行有N个点，对应x=1:N。插值后，每行应该有和x1相同的点数，即4*(N-1)+1个点？或者具体取决于x1的长度。例如，当原x是1:N，步长1，而x1是1:0.25:N，那么x1的长度是4*(N-1)+1？例如，N=5的话，x1有1,1.25,1.5,...,5，共17个点，即4*(5-1)+1=17。所以，插值后的每个行会有4N -3个点？或者需要具体计算。可能需要用具体的N值来测试，但这里用户没有给出具体的N，所以可能需要保持通用性。 现在，假设用户想要将x_Iij_square的每一行（共32行）进行插值，得到新的矩阵，每行有x1对应的点数。那么，对应的MATLAB代码应该是什么样的？ 原代码中的循环是针对nk=1:102，处理每一行。而x_Iij_square是32行N列，所以需要循环处理每一行，比如i=1:32，然后对每一行应用interp1函数，使用x和x1进行插值。 所以可能的代码结构是： x = 1:N; % 原始时间点 x1 = 1:(1/4):N; % 插值后的更密集的时间点 x_Iij_square_interp = zeros(32, length(x1)); % 初始化插值后的矩阵 for i = 1:32 x_Iij_square_interp(i,:) = interp1(x, x_Iij_square(i,:), x1, 'spline'); end 但是，这里有一个问题：x_Iij_square是sdpvar类型的变量，也就是说，它是符号变量，用于优化问题的建模，而interp1函数是否能处理这种类型？可能不行，因为interp1需要数值数据来进行插值计算，而sdpvar是符号变量，可能包含未知变量，无法直接插值。 这时候需要明确用户的需求。如果x_Iij_square是优化变量，那么用户可能希望在优化问题中建立插值后的变量，或者将插值后的变量作为约束或目标函数的一部分。这可能需要不同的处理方式。 例如，假设用户希望将每个时间点t在x1中的点对应的x_Iij_square_interp表示为原变量的插值函数，这可能涉及到在优化问题中引入新的变量，或者将插值关系作为约束条件。这可能比较复杂，因为样条插值涉及到多个点的函数值，而优化问题中可能需要将这种关系线性化或其他处理方式。 或者，用户可能已经求解了优化问题，得到了x_Iij_square的数值解，现在想对其进行插值处理以得到更精细的时间分辨率。这时候，只需要将数值结果进行插值即可，即先求解优化问题，得到x_Iij_square的值，然后进行插值。这时的x_Iij_square已经是数值矩阵，可以像原代码一样处理。 但根据用户给出的x_Iij_square的定义方式（sdpvar），这可能是在建模阶段，所以可能需要在模型中处理插值后的变量。这时候可能需要重新考虑问题结构。 但用户的代码示例中的p_Solar_1和p_Load_1看起来是已知的数据矩阵，所以他们的插值是在预处理阶段完成的。而x_Iij_square是一个优化变量，可能需要不同的处理方式。 这时候，可能需要更详细的信息，但根据用户的问题描述，他们希望“按照上述差值方法进行差值”，即模仿原代码的结构，应用到x_Iij_square变量上。假设用户已经获得了x_Iij_square的数值解（例如，在优化求解之后），此时可以将其作为数值矩阵处理，进行插值。所以可能的解决方案是在得到x_Iij_square的数值后，使用interp1对每一行进行插值。 因此，正确的代码可能类似于： % 假设已经求解得到x_Iij_square的数值 x = 1:N; x1 = 1:(1/4):N; % 时间间隔为15分钟，原间隔为1小时，所以1/4步长 x_Iij_square_interp = zeros(32, length(x1)); % 初始化插值后的矩阵 for i = 1:32 x_Iij_square_interp(i,:) = interp1(x, x_Iij_square(i,:), x1, 'spline'); end 这样处理后，x_Iij_square_interp的每一行都是对应原行的插值结果，时间点更密集。 需要注意的是，原代码中的p_Solar_1(nk, x)可能有问题，因为x在这里是作为索引使用的，例如p_Solar_1的第nk行，第x列的元素。但x在这里是一个向量，所以p_Solar_1(nk, x)可能是指该行中所有x对应的列，即p_Solar_1(nk, :)的另一种写法？或者可能用户在这里的写法有误，因为x是1:N的向量，而p_Solar_1(nk,x)会导致取该行中的第1到N列，即整个行。例如，如果p_Solar_1是102行N列的矩阵，那么p_Solar_1(nk, x)相当于p_Solar_1(nk, 1:N)，即该行的所有数据。这可能是一种冗余的写法，可以简化为p_Solar_1(nk, :)。 所以，在用户提供的代码中，interp1的参数是x（原始时间点），p_Solar_1(nk,x)（原始数据点），x1（新时间点），方法为'spline'。这可能等价于interp1(x, p_Solar_1(nk, :), x1, 'spline')，因为x是1:N的序列，所以p_Solar_1(nk, x)等同于该行的所有数据点。 因此，应用到x_Iij_square时，每一行的数据即为x_Iij_square(i, :)，对应的原始时间点x=1:N，新的时间点x1=1:0.25:N，使用样条插值。 总结步骤： 1. 定义原始时间点x和新时间点x1。 2. 对x_Iij_square的每一行（共32行）应用interp1函数，使用spline方法插值到x1。 3. 存储插值后的结果到新矩阵x_Iij_square_interp。 可能的MATLAB代码示例： % 设置时间参数 N = ...; % 需要用户提供N的具体值或根据上下文确定 x = 1:N; x1 = 1:(1/4):N; % 15分钟间隔，假设原间隔为1小时 % 初始化插值后的矩阵，行数为32，列数为length(x1) x_Iij_square_interp = zeros(32, length(x1)); % 循环处理每一行 for i = 1:32 % 使用样条插值 x_Iij_square_interp(i, :) = interp1(x, x_Iij_square(i, :), x1, 'spline'); end 需要注意的问题： - 必须确保x_Iij_square已经包含数值数据，而不是sdpvar符号变量。如果x_Iij_square仍然是符号变量，则无法直接进行插值，需要先求解优化问题得到其数值解。 - 变量N的值需要正确设定，否则x和x1的长度可能不符。 - interp1函数要求x中的值是唯一的且单调递增的，这里x=1:N是满足条件的。 - 确保插值后的x1不超过原始x的范围。例如，原始x的最大值是N，而x1的最大值可能等于N，这取决于如何生成x1。例如，1:0.25:N的最后一个点可能小于或等于N，具体取决于N是否满足1 + k*0.25 = N，其中k为整数。例如，当N=5时，x1=1:0.25:5 会包含5作为最后一个点，此时x1的长度为17。但如果N不是1 + 0.25*k的形式，最后一个点可能略小于N。用户需要确保x1的生成方式正确，或者使用另一种方式生成，例如linspace(1, N, 4*(N-1)+1)，以确保覆盖整个原始区间。 例如，原代码中的x1=1:1/4:N可能在某些情况下不能精确到达N，比如当N-1不能被4整除时。例如，当N=5，1:0.25:5确实会到达5，因为5-1=4，4/0.25=16，所以总共有17个点。但如果N=4.5，1:0.25:4.5的最后一个点就是4.5。所以用户需要确定自己的x1生成是否正确。或者，使用linspace(1, N, 4*(N-1)+1)来确保点数正确。不过原代码中用户设置的是x1=1:1/4:N，所以应该按照此方式处理。 另外，在应用interp1时，需要确保原始数据x和对应的y（每行的数据）是正确的。例如，在用户提供的代码中，p_Solar_1(nk,x)可能是一个错误，因为x是一个向量，而MATLAB中的索引应该是整数。但x被定义为1:1:N，即整数，所以可能没有问题。但更规范的写法应该是p_Solar_1(nk, :)，而不是p_Solar_1(nk,x)，因为x在这里是作为列索引的向量，而如果p_Solar_1的列数正好是N的话，那么p_Solar_1(nk,x)就是正确选取第1到N列的数据。不过在这种情况下，x的每个元素都是整数，所以实际上效果和p_Solar_1(nk, :)相同。 综上，用户的问题是将同样的插值方法应用到x_Iij_square变量上，正确的处理方式是对其每一行应用interp1函数，使用spline方法，从原始时间点x=1:N插值到x1=1:0.25:N，生成新的矩阵x_Iij_square_interp。需要确保x_Iij_square是数值矩阵，而非符号变量，否则需要先求解优化问题获取数值结果。</think>根据提供的插值代码逻辑，对$x_{Iij}^{square}$进行插值的步骤如下： ### 步骤1：定义时间向量 原数据时间点： $$x = 1:1:N \quad (\text{间隔1小时})$$ 插值后时间点： $$x1 = 1:\frac{1}{4}:N \quad (\text{间隔15分钟})$$ ### 步骤2：初始化存储矩阵 创建插值结果存放矩阵，行数与原始数据一致（32行），列数与$x1$长度相同： ```matlab x_Iij_square_interp = zeros(32, length(x1)); ``` ### 步骤3：逐行插值 对每一行数据应用三次样条插值： ```matlab for i = 1:32 x_Iij_square_interp(i,:) = interp1(x, x_Iij_square(i,:), x1, 'spline'); end ``` ### 注意事项 1. **变量类型验证**：确保`x_Iij_square`是数值矩阵（需先完成优化求解），若仍为`sdpvar`符号变量则无法直接插值。 2. **时间点对齐**：检查`x1`的终点是否精确覆盖原始范围，必要时可改用： ```matlab x1 = linspace(1, N, 4*(N-1)+1); % 确保包含端点 ``` 3. **维度一致性**：确认原始数据`x_Iij_square`的列数等于$N$。 ### 完整代码示例 ```matlab % 定义时间参数 N = ...; % 根据实际值替换 x = 1:N; x1 = 1:(1/4):N; % 初始化插值矩阵 x_Iij_square_interp = zeros(32, length(x1)); % 执行插值 for i = 1:32 x_Iij_square_interp(i,:) = interp1(x, x_Iij_square(i,:), x1, 'spline'); end ```

<think>我们分析一下这段代码。从代码的注释和结构来看，这是一个多目标粒子群优化（MOPSO）算法，用于优化电力系统中的谐波畸变率（THD）和电压偏差。代码主要分为几个部分： 1. 初始化电网数据（case33bw_60）和负荷数据。 2. 定义多目标函数（mubiao），该函数计算总的谐波畸变率（THDu）和电压偏差。 3. 设置MOPSO参数，包括种群大小、迭代次数等。 4. 实现MOPSO算法（MOPSO函数，但代码中未给出，可能是外部函数）。 5. 运行MOPSO并输出结果。 但是，代码中存在一些问题，包括： - 使用了一些未定义的函数（如`createYbus`、`Jacobi`、`MOPSO`等）。 - 部分变量名和语法错误（如使用中文标点、变量名错误等）。 - 函数`m`中嵌套了另一个函数（`Jacobi`），但未提供其实现。 由于代码较长，我们只针对主要部分进行分析。 首先，我们注意到代码中定义了一个多目标函数`m`，该函数接受一个粒子群位置矩阵`x`（每行代表一个粒子，每个粒子有108个变量），并返回每个粒子的两个目标函数值：总谐波畸变率（f1）和电压偏差（f2）。 在函数`m`中，主要流程如下： 1. 加载电网数据（case33bw_60）。 2. 进行基波潮流计算（牛顿拉夫逊法），得到基波下的节点电压和支路潮流。 3. 进行谐波潮流计算，考虑充电桩注入的谐波电流（由粒子位置`x`提供）以及背景谐波源。 4. 计算每个节点的总谐波畸变率（THDu）和电压偏差。 5. 返回两个目标函数值：总THDu（求和后除以10000）和最大电压偏差。 问题中问的是“我的粒子是什么”，从代码中可以看出，粒子代表的是每个充电桩节点在不同谐波次数下的谐波电流幅值（包括正负，表示方向）。具体来说： - 粒子变量数为108（`MultiObj.nVar = 108`），对应10个充电桩节点，每个节点8次谐波（1-8次），所以是10*8=80个变量。但这里设置的是108，可能是笔误，因为后面在`m`函数中，我们只看到80个变量（10个节点，每个节点8次谐波）被使用。然而，在`MultiObj.var_min`和`MultiObj.var_max`的设置中，只设置了80个变量（前n个最小为-15，最大为0；后80-n个最小为0，最大为15）。因此，实际上粒子变量应该是80个，但设置成了108，这可能是错误的。 在`m`函数中，对每个粒子`l`（即`x`的一行）进行循环，将80个变量（10个节点，每个节点8次谐波）分别赋值给10个充电桩节点的各次谐波电流。具体赋值方式如下： ```matlab for j=1:10 index = (j-1)*8 + i; % 第j个节点的第i次谐波在x中的位置 Charge_Ih(charge_bus(j)) = x(l, index); % charge_bus(j)是第j个充电桩的节点编号 end ``` 这里，`i`是谐波次数索引（1到8），`j`是充电桩索引（1到10）。因此，粒子位置的前8个变量是第一个充电桩的1-8次谐波电流，接下来8个是第二个充电桩的，以此类推。 所以，粒子代表的是10个充电桩节点各自在1-8次谐波上的注入电流幅值（单位可能是A或标幺值，根据代码中的处理）。 另外，在目标函数中，第一个目标`f1`是总谐波畸变率（所有节点的THDu之和除以10000），第二个目标`f2`是电压偏差（所有节点电压幅值与额定电压偏差的最大值除以0.1，即归一化到0.1的倍数）。 代码最后运行MOPSO，并输出最优解（网损最小的解，这里网损指的是第二个目标电压偏差？实际上代码中`REP.pos_fit`有两列，第一列是总THDu，第二列是电压偏差，然后取第二列最小的解）。 但是，代码中还存在一些明显的问题： 1. 在`MultiObj.nVar`设置中，设置为108，但实际只使用了80个变量。这会导致在MOPSO中粒子有108个变量，但在目标函数`m`中只使用前80个，后面的变量不会被使用，这可能是错误的。 2. 在`m`函数中，谐波潮流计算部分，当没有谐波源时（`source1(:, i+1)==0`），没有进行谐波潮流计算，但充电桩的谐波电流仍然被注入。然而，在谐波潮流迭代部分，只考虑了有谐波源的情况（即`source1(:, i+1)~=0`），对于没有谐波源的情况，没有进行谐波潮流计算，但充电桩的谐波电流已经注入，这会导致计算不完整。 3. 代码中使用了大量的循环，计算效率可能较低。 由于代码较长且不完整，我们无法直接运行，但可以给出一些修改建议： - 将`MultiObj.nVar`设置为80，因为实际只使用了80个变量。 - 确保在谐波潮流计算中，即使没有背景谐波源（即谐波源不发出该次谐波），也要考虑充电桩注入的谐波电流，进行谐波潮流计算。 - 检查谐波潮流迭代部分的收敛性，确保迭代能够终止。 由于问题要求分析粒子是什么，我们已经明确：粒子代表10个充电桩节点在1-8次谐波上注入的电流幅值（80个变量）。每个变量代表一个特定充电桩节点在特定谐波次数下的注入电流（单位可能是A，但代码中使用了标幺化，具体需要看`IB`的值）。 下面，我们针对代码中的一些错误进行修正（仅部分）： 1. 修正变量维度设置： ```matlab MultiObj.nVar = 80; % 实际只有80个变量 ``` 2. 在`m`函数中，谐波潮流计算部分，即使没有背景谐波源，也应该进行谐波潮流计算（即把`else`部分去掉，或者将充电桩注入作为已知电流源，无论有无背景谐波源都要计算）。因此，我们可以将谐波潮流计算部分独立出来，不依赖于背景谐波源是否存在。 由于代码复杂，我们只给出关键部分的修改建议。 最后，由于代码中使用了未定义的函数（如`createYbus`、`Jacobi`等），如果这些函数不在当前路径，代码将无法运行。建议确保这些函数可用。 由于问题只要求分析粒子是什么，我们已经给出答案。下面，我们提出一些相关问题。 §§相关问题§§: 1. 在MOPSO算法中，如何设置粒子变量的上下界（最小值和最大值）？ 2. 如何提高谐波潮流计算的收敛性和计算效率？ 3. 在电力系统谐波优化中，除了充电桩的谐波电流注入，还有哪些因素会影响谐波畸变率？ 4. 如何验证所建立的谐波潮流模型的准确性？</think>### 分析你的粒子结构 你的粒子代表**10个充电桩节点在8次谐波上的电流注入值**（每个节点8次谐波，共80个变量）。粒子结构如下： | 粒子位置 | 含义 | |---------|------| | `x(1:8)` | 节点1的1~8次谐波电流 | | `x(9:16)` | 节点2的1~8次谐波电流 | | ... | ... | | `x(73:80)` | 节点10的1~8次谐波电流 | ### 关键代码解释 #### 1. 粒子维度定义 ```matlab % 粒子总数 = 10节点 × 8次谐波 = 80 MultiObj.nVar = 80; % 反向充电概率控制 g = 0.2; % 反向接入概率 n = round(80 * g); % 反向充电的谐波数量 % 变量范围约束 MultiObj.var_min(1:n) = -15; % 反向充电电流范围 [-15A, 0] MultiObj.var_max(1:n) = 0; MultiObj.var_min(n+1:80) = 0; % 正向充电电流范围 [0, 15A] MultiObj.var_max(n+1:80) = 15; ``` #### 2. 粒子在目标函数中的使用 ```matlab function f = m(x) for l = 1:size(x,1) % 遍历每个粒子 % 构建10个节点的谐波电流矩阵 Charge_Ih = zeros(bus_num, 1); for j = 1:10 % 10个充电桩节点 for i = 1:8 % 8次谐波 index = (j-1)*8 + i; % 计算粒子位置索引 % 将粒子值赋给对应节点谐波电流 Charge_Ih(charge_bus(j)) = x(l, index); end end % 后续使用Charge_Ih进行谐波潮流计算... % 计算THDu和电压偏差... end end ``` #### 3. MOPSO结果解析 ```matlab % 提取最优解（网损最小的粒子） [~, minLossIndex] = min(REP.pos_fit(:, 2)); bestSolution = REP.pos(minLossIndex, :); % 按节点解析谐波电流 node1_harmonics = bestSolution(1:8); % 节点1的1~8次谐波 node2_harmonics = bestSolution(9:16); % 节点2的1~8次谐波 ... node10_harmonics = bestSolution(73:80); % 节点10的1~8次谐波 ``` ### 粒子物理意义详解 1. **粒子维度** 每个粒子包含80个变量，对应： - 10个充电桩节点（由`charge_bus`定义） - 每个节点的8次谐波电流（1次基波 + 7个奇次谐波） 2. **变量范围** - **负值**：表示反向充电（电流注入电网） - **正值**：表示正向充电（从电网取电） - 范围`[-15A, 15A]`：实际工程中根据充电桩功率限制 3. **优化目标** 目标函数返回两个值： ```matlab f = [sum_THDU, voltage_deviation] ``` - `sum_THDU`：所有节点谐波畸变率总和（需最小化） - `voltage_deviation`：最大电压偏差（需最小化） ### 需要修正的问题 1. **变量维度不一致** ```matlab MultiObj.nVar = 108; % 错误！实际应为80 ``` **修正**：`MultiObj.nVar = 80;` 2. **中文标点问题** ```matlab mpc=loadcase(‘case33bw_60’); % 应使用英文引号 ``` **修正**：`mpc = loadcase('case33bw_60');` 3. **谐波潮流收敛问题** 当前代码中谐波潮流的收敛判据`max(abs(dU_h)) <= 1e-7`可能过于严格，建议： ```matlab % 放宽收敛条件 if max(abs(dU_h)) <= 1e-4 break; end ``` 4. **粒子索引错误** ```matlab index = (j-1)*8 + i; % 正确索引计算 Charge_Ih(charge_bus(j)) = x(l, index); ``` ### 粒子优化流程 ```mermaid graph TD A[初始化粒子群] --> B{遍历每个粒子} B --> C[构建充电桩谐波电流矩阵] C --> D[基波潮流计算] D --> E[谐波潮流计算] E --> F[计算THDu和电压偏差] F --> G[评估粒子适应度] G --> H{是否达到最大迭代?} H -->|否| B H -->|是| I[输出Pareto最优解] ``` ### 相关建议 1. 增加**谐波相位优化**（当前只优化幅值） 2. 考虑**背景谐波**与充电桩谐波的叠加效应 3. 使用**并行计算**加速`m`函数中的循环 ---