预处理计算方法
ILU预处理Newton—Krylov方法的潮流计算
ILU预处理Newton—Krylov方法的潮流计算作者：廖小兵王文超李奔来源：《计算技术与自动化》2015年第04期摘要：由于电力系统修正方程组具有高维、稀疏的特点，本文提出将预处理Krylov子空间方法应用于潮流修正方程组的求解，形成预处理NewtonKrylov的潮流计算方法。
结合ILU 预处理方法，比较了最常用的3类NewtonKrylov方法求解潮流方程的计算效果。
通过对IEEE30、IEEE118、IEEE300 和3个Poland大规模电力系统进行潮流计算，结果表明：3类NewtonKrylov方法是电力系统潮流计算的有效方法，呈现出良好的收敛特性和计算效率。
关键词：潮流计算；修正方程组；ILU预处理；NewtonKrylov方法中图分类号：TM744 文献标识码：A1引言潮流计算是电力系统分析中最古老的经典课题之一。
传统的电力系统潮流计算通常以牛顿法为主[1]。
牛顿法是求解非线性代数方程组的有效方法之一，它将非线性方程组的求解转化为线性代数方程组的求解，但由于每次迭代后都需更新雅可比矩阵的元素，导致每次都需求解高维的线性代数方程组。
传统的直接法，如Gauss消去法，LU分解等，计算量和存储量较大，且固有的前推回代过程难以并行[2-3]。
迄今为止，越来越多的国内外研究人员在电力系统潮流计算中采用NewtonKrylov方法求解潮流方程[4-8]。
NewtonKrylov方法是在不精确牛顿法的基础上，结合Krylov子空间迭代法，形成的一类新的求解非线性方程组的数值方法。
这类方法结合了Newton方法的良好收敛特性，以及Krylov子空间方法的存储量少、计算量小、易于并行等优点[9]，非常适合并行求解大规模的非线性方程组问题[10]。
文献[4]首次将Krylov子空间法中的GMRES方法应用于潮流计算中。
文献[5]将此类迭代法与不精确牛顿法相结合（NewtonGMRES），同时采用不同的预处理方法，对两个大规模电力系统进行了对比分析计算。
基于预处理的非稳态流场计算方法
中圈分类 号 ：４ ０ Ｖ ３ 文献标识码 ： Ａ 文章 编号 ：０６２ ９ （ｏ６ ０ ｌ９ ｌｏ ＿７ ｆ ｓｏ ａｔ ｓ Ｂ ｒｎ ｎｖ ｆ ｅｏａｔ ｓａｄＡｔ ｎｕｉ ， ｅ ｉｇ １０ ８ ，ｈｎ ） Ｓ ｈｏ ｏ ｔ ｎｕｉ ， ｅ ｉｇＵ ｉ．ｏ ｒｎｕ ｃ ｎ ｓ ａｔ ｓ Ｂ ｒｎ ０ ０ ３ Ｃ ｉ Ａ ｒ ｃ Ａ ｉ ｏ ｒ ｃ ａ
Ａｂ ｔａ ｔ Ｔ ｅａ ｃ ｒｃ ｎ ｏ ｖ ｒｅ ｃ ｈ ｒｃｅｉｔ ｓｏ Ｓ ｏ ｕｅ ｑ ａｉｎ ｓｄ ｆｒｃ ｃ ｌｔ ｇｌｗ Ｍａｈ ｎ ｍｂ ｒｆ ｗ ｓｒ ｃ ： ｈ ｃ ｕａ ｙａ ｄ ｃ ｎ ｅｇ ｎ ｅｃ ａａｔｒ ｉ ｆＮ— ｒＥ ｒｅ ｕ ｔ ｓｕ ｅ ｏ ａ ｕａ ｎ ｏ ｃ ｕ ｅ ｏ ｓｃ ｌ ｏ ｌ ｉ ｌ ｉ ｄ ｗ ｒ ｉ ｒ ｅ ｙｕ ｉｇｐ ｅｏ ｄｔ ｎｎ ｔｏ ．Ｔｍｅｐｅ ｉｉ ｆ ｌ ｕａｉ ｕｔ ｏ ｎ ｔ ｄ ｏ ｅｄＷ Ｓｅ ｓｒ ｙ ｕ ｆ ｌ ｅ ｍｐ ｖ ｄ ｂ ｓｎ ｒｃ ｎ ｉｏ ｉｇｍｅｈ ｄ ｉ ｒｃｓ ｎｏ ａｃ ｌｔ ｎ ｒｓ ｌ ｆｕ ｓ ａ ｙｆ ｗ ｆ ｌ ａ ｎ ｕ ｄｂ － ｅ ｅ ｏ ｉ ｏ ｃ ｏ ｅ ｓ ｅ ｌ ｉ ｅ
ｔ ｃ ｅ ｄｃ ｎ ｅ ｔ ｎｆ ｘｕ ｅ ｉ ｓ ｈ ｍｅａ ｏ ｖ ｃｉ ｕ ｓ ｄＡＵＳ ＋ ａ ｄＭ Ｕ Ｃ ｃ ｅ Ｓｓｌｅ ｔ ｖ ｒ ｈｒｃ ｌ ｍｅｓｅ ．Ｔ ｒ￣， ａｃ ｌｔ ｎｏ ｍｅ ｎ ｏ ｌ Ｍ ｎ Ｓ Ｌ ｓｈ ｍｅＷａ ｏｖ ｄａ ｅｙ ｐ ｙ ｉａ ｅ ｓ ｉ ｔ ｔｐ ｈ ｏ ｈｃ ｕａ ｏ ｆ ｌ ｉ ｆｒ ａｄ ｆｃｎ ｔｐ ｆ ｗ ａ ｄｃｒｕａ ｙｉｄ ｒ ｈ ｌｏｉ ｍ Ｓｐ ｖ ｎ ｔ ｅ ｏ ｉｅａ ｐ ｉａ ｉｔ ｎ ｅｉ ｉｔ． ｏｗ ｒ — ｉｇｓｅ ｏ ｉ ｌｒ ｃ ｌ ｅ ，ｔｅａ ｒｈ Ｗａ ｒ ｅ ｏｂ ｆｗ ｄ ｐ ｌ ｂｌｙａ ｄｒｌ ｌ ａ ｌ ｎ ｃ ｎ ｇ ｔ ｏ ｃ ｉ ｂ ｙ ａ ｉ Ｋｅ ｒ ：ｕ ｌｔ ｔｐ；ｐ ｃ ｎ ｉｏ ｉｇ ； ｎｔａ ｙｆ ｗ ｙ ｗｏｄｓ ｄ ａ－ｉ ｓｅ ｍｅ ｅ ｒ ｏ ｄｔ ｎｎ ｕ ｓｅ ｄ ｏ ｉ ｌ
小水电生态径流计算数据预处理方法
在置信 区间 ９ ５ ％以内的合理的数 据 ，三角形点表
示 波动 异常 的数据 。
做 出 回归 曲线之 后 ，通过描 绘残 差 图来进行 回
归诊 断 。所 描 绘 的 残 差 图是 以残 差 ｒ 或 ｅ为 纵 坐
标 ，以任 何其他 的量 为横 坐标 的散点 图。通过残 差 图的分析 ，一方 面 可 以检 验求解 出的 回归 曲线是 否 符合 实 际 的分 布 ，另一方 面可 以检 验数据 点是 否落 在 执行 区 间之 内 ，若为 异常 点 （ 历 史水 文数 据 中误 差 较大 的水文 记 录 ）或 者强影 响点 （ 限水 文 年限 中
小 ，上述 的两类 数 据将 可能 导致生 态灾 难 ，因此 在 小水 电生 态径 流 的计 算 中对 于历史 水 文数据 的预 于此 ，提 出基于线 性 回归统 计对 历史 水 文数 据做 预处 理 ，筛 选剔 除 波动 幅度异 常 的数据 并用 线
性回归值替代 ，来提高最后的生态径流计算数值 的
图 １ ） 。
数 量级 较 大 的天然径 流 ，历 史水 文数 据 的可 能误差 以及极 限水 文 年 的极 限水文 数据 不会 直接 影 响到 生
态径流最后 的计算值 。但对于径流量数量级较小的
天然 径 流 ，尤 其 是 小 水 电 所 处 的径 流 流 量 普 遍 较
收 稿 日期 ：２ ０ １ ３ —０ ７ —２ ０
小水电 ２ ０ １ ３ 年第 ５ 期 （ 总第 １ ７ ３ 期）
农村 水 电增 效扩容 改造
小水 电生态径 流计 算数 据预 处理 方法
周 李 智 王万 良 徐新 黎 赵 燕伟 （ 浙 江工业 大 学 浙 江杭 州 ３ １ ０ ０ ２ ３ ）
python数据预处理的方法
python数据预处理的方法
Python中常用的数据预处理方法有以下几种：
1. 数据清洗：对数据进行去重、缺失值处理、异常值处理、数据类型转换等操作，以保证数据的准确性和完整性。
2. 特征选择：根据特征重要性和相关性等指标，选择最具代表性和区分度的特征，以提高模型的准确性和泛化能力。
3. 特征缩放：对特征进行缩放操作，使得不同特征之间具有相同的尺度，以便于模型的训练和优化。
4. 特征编码：将文本型特征进行编码，以便于模型的处理和计算。
常见的编码方法包括one-hot编码、标签编码、二进制编码等。
5. 数据归一化：将数据进行归一化处理，以使得数据分布在[0,1]或[-1,1]之间，以便于模型的训练和优化。
6. 数据降维：对高维数据进行降维处理，以减少特征数量和计算复杂度，常用的降维方法包括主成分分析法（PCA）、线性判别分析法（LDA）等。
7. 数据集划分：将数据集划分为训练集、验证集和测试集，以便于模型的训练、评估和优化。
以上是Python中常用的数据预处理方法，不同的数据预处理方法可以根据具体情况进行选择和组合。
带状矩阵方程预处理方法
带状矩阵方程预处理方法
嘿，咱今儿就来唠唠带状矩阵方程预处理方法。
这玩意儿啊，就好
像是一把解开复杂难题的钥匙。
你想想看，矩阵就像是一个大迷宫，而带状矩阵呢，就是这个迷宫
里有规律可循的那部分。
处理带状矩阵方程，就像是在迷宫中找到一
条最快捷的路。
预处理方法呢，就好比给这条路先清理清理障碍，让我们能更顺利
地走过去。
它可以让计算变得更高效，就像给汽车加上了超强马力。
比如说，我们在处理一些大型数据的时候，要是没有好的预处理方法，那可就像没头苍蝇一样乱撞啦。
但有了它，就能快速准确地找到
答案。
咱再打个比方，这预处理方法就像是给运动员做的热身准备，只有
充分热身了，才能在比赛中发挥出最好的水平呀。
它能让整个计算过
程变得更流畅，更轻松。
那怎么选择合适的预处理方法呢？这可得好好琢磨琢磨。
不同的问
题可能需要不同的方法，就像不同的病症需要不同的药一样。
有时候
可能需要试试几种方法，才能找到最适合的那个。
而且啊，这可不是一劳永逸的事儿。
就像人的习惯会变一样，问题
也可能会变，那预处理方法也得跟着变呀。
不能死脑筋，得灵活点儿。
有时候可能会遇到一些特别难搞的情况，就像爬山遇到了陡峭的山峰。
但别怕呀，咱有预处理方法这根“登山杖”呢，总能找到办法爬上去的。
总之呢，带状矩阵方程预处理方法是个很重要的东西。
它能让我们在解决复杂问题的时候更得心应手，就像有了一个得力的助手。
咱可得好好研究研究它，让它为我们的学习和工作发挥最大的作用呀！可别小瞧了它哦！。
大规模数据分析中的计算方法
大规模数据分析中的计算方法在当今数字化时代，数据分析已经成为商业、科学甚至社会决策的重要基础。
然而，大规模数据的处理和分析需要高效的计算方法来实现。
本文将探讨现代大规模数据分析中常用的计算方法，包括数据预处理、机器学习和深度学习等方面。
一、数据预处理的计算方法在大规模数据分析中，数据预处理是一个必不可少的步骤，因为数据的质量和准确性直接影响到后续分析的可靠性。
数据预处理的主要目的是去除重复数据、异常数据，填补缺失值，并对数据进行标准化和归一化等操作。
为了实现数据预处理，常用的计算方法包括以下几个方面：1. 去重。
可以使用哈希表或者位图等方法实现快速查找与去重。
以哈希表为例，在将数据插入哈希表时，如果发现哈希表中已经存在该数据，则删除该数据。
2. 数据填充。
数据填充可以使用线性插值、KNN填充等方法。
以KNN填充为例，首先选择一个KNN算法，然后选择适当的K值，计算出缺失值周围的K个邻居，最后使用这K个邻居的平均值或中位数来填充缺失值。
3. 标准化和归一化。
标准化和归一化的目的是将数据转换为一定范围内的值，以便更好地进行分析。
常用的方法包括Min-Max Scaling以及Z-Score Scaling。
其中，Min-Max Scaling将数据转换到[0,1]之间，而Z-Score Scaling将数据转换为均值为0，方差为1的正态分布。
二、机器学习中的计算方法机器学习是一种通过训练从数据中提取规律的方法，是大规模数据分析中最常用的方法之一。
常见的机器学习算法包括决策树、支持向量机、朴素贝叶斯、KNN等。
这些算法通过学习输入数据和答案之间的关系，并在更大的数据集上进行测试和预测。
对于这些算法，计算方法包括以下几个方面：1. 特征提取。
特征提取是机器学习过程中最关键的步骤之一。
它涉及到选择合适的特征，将复杂的输入数据转换为简单的数值向量。
例如，在进行图像识别时，特征提取就可以将每个图像转换为一组数字，表示其纹理、形状和颜色等特征。
预处理方法
第五讲预处理方法1.预处理Krylov子空间方法2.基于矩阵分裂的预处理技术3.基于不完全分解的预处理子4.基于稀疏近似逆的预处理子5.其他预处理技术为什么预处理Nothing will be more central to computational science in the next century thanthe art of transforming a problem that appears intractable into another whosesolution can be approximated rapidly.For Krylov subspace matrix iterations,this is preconditioning.—Trefethen&Bau,1997.•为什么预处理：改善迭代法的收敛性与可靠性.•预处理的根本目的：将难求解的原问题转化成等价的容易求解的新问题.•1968年,Evans首先提出“preconditioning”,并考虑了不完全LU分解预处理.•预处理技术的研究是目前科学计算领域中的重要研究课题之一.1预处理Krylov子空间方法1.1怎么做预处理1.2预处理CG方法1.3预处理GMRES方法1.4左、右预处理方式的最优性质1.5Flexible GMRES方法1.6预处理BiCGSTAB方法1.1怎么做预处理预处理方法考虑线性方程组Ax=b,A∈R n×n非奇异,b∈R n.(5.1)在方程组的两边同时左乘一个非奇异矩阵P的逆,则可得P−1Ax=P−1b(5.2)这个新方程组就是预处理后的方程组,P就称为预处理子(preconditioner).当Krylov子空间方法被用于求解(5.2)时,就称为预处理Krylov子空间方法.好的预处理子理论上讲,任何一个非奇异矩阵都可以作为预处理子.但一个好的预处理子P通常需满足下面两个要求:(1)P容易构造且P z=r容易求解→实用(2)大大提高迭代法的收敛速度→有效（P−1A具有更小的条件数,更好的特征值分布,次数更低的最小多项式）三种不同的预处理方式预处理方式(5.2)称为左预处理,即从左边做预处理.相应地,可以对系数矩阵右乘P−1,这就是右预处理,即AP−1u=b,x=P−1u(5.3)另外,如果P可以分解成两个矩阵的乘积,或者P是由两个矩阵的乘积构成,即P=P L P R,则可以用下面的方式对原方程组进行预处理P−1 L AP−1Ru=L−1b,x=P−1Ru(5.4)这就是两边预处理三种预处理方式的几点说明以上三种方式预处理后的系数矩阵分别为P−1A,AP−1和P−1L AP−1R:•具有相同的特征值分布(三个矩阵是相似的)•如果A和P都是对称正定的,则使用共轭梯度法求解时,这三种方式的预处理效果基本上是一样的.•对于非对称(特别是非正规)情形,效果可能会相差很大.三种预处理方式的几点说明（续）1.2预处理CG方法PCG设A对称正定,为了保证预处理后的系数矩阵仍然是对称正定的,通常要求预处理子P也是对称正定的.考虑使用两边预处理方式.设P的Cholesky分解为P=LL⊺.于是我们得到下面的预处理方程组L−1A(L⊺)−1u=L−1b,x=(L⊺)−1u.(5.5)近似解记为u (k ),预处理残量记为˜r k ≜L −1b −L −1A (L ⊺)−1u (k ).算法1.1两边预处理CG 方法1:给定初值x (0),计算r 0=b −Ax (0)2:令˜r 0=L −1r 0,˜p 1=˜r 03:for k =1,2,...do 4:ξk =(˜r k −1,˜r k −1)(L −1A (L ⊺)−1˜p k ,˜p k )5:u (k )=u (k −1)+ξk ˜p k 6:˜r k =˜r k −1−ξk L −1A (L ⊺)−1˜p k 7:µk =(˜r k ,˜r k )(˜r k −1,˜r k −1)8:˜p k +1=˜r k +µk ˜p k 9:end for原方程组残量PCG得到的是预处理后的近似解u(k)和残量˜r k,因此我们还需要考虑原方程组的近似解和残量的计算.对上述算法中进行适当的改写.首先引入一个辅助变量p k:p k≜(L⊺)−1˜p k,k=1,2,....于是有p k+1=(L⊺)−1p k+1=(L⊺)−1˜r k+µk(L⊺)−1˜p k=(L⊺)−1˜r k+µk p k.由˜r k的表达式可知,原方程组的残量r k=L˜r k,所以有r k=L˜r k=L(˜r k−1−ξk L−1A(L⊺)−1˜p k)=r k−1−ξk Ap k,p k+1=(L⊺)−1˜r k+µk p k=(L⊺)−1L−1r k+µk p k=P−1r k+µk p k,x(k)=(L⊺)−1u(k)=(L⊺)−1(u(k−1)+ξk˜p k)=x(k−1)+ξk p k,其中ξk=(L−1r k−1,L−1r k−1)(L−1A(L⊺)−1˜p k,˜p k)=(r k−1,(L⊺)−1L−1r k−1)(A(L⊺)−1˜p k,(L⊺)−1˜p k)=(r k−1,P−1r k−1)(Ap k,p k),µk=(L−1r k,L−1r k)(L−1r k−1,L−1r k−1)=(r k,P−1r k)(r k−1,P−1r k−1).记z k≜P−1r k,则可得下面的预处理CG方法.算法1.2预处理CG方法(PCG)1:给定初值x(0),(相对)精度要求ε>0和最大迭代步数IterMax 2:计算r0=b−Ax(0)和β=∥r0∥23:令z0=P−1r0,p1=z04:计算ρ=r⊺0z05:for k=1to IterMax do6:q k=Ap k7:ξk=ρ/(p⊺k q k)8:x(k)=x(k−1)+ξk p k9:r k=r k−1−ξk q k10:relres=∥r k∥2/β11:if relres<εthen12:break13:end if14:ρ0=ρ15:z k=P−1r k16:ρ=r⊺k z k17:µk=ρ/ρ018:p k+1=z k+µk p k19:end for20:if relres<εthen21:输出近似解x(k)及相关信息22:else23:输出算法失败信息24:end if左预处理方式事实上,PCG也可以通过左预处理方式导出.考虑左预处理方程组P−1Ax=P−1b.(5.6)易知P−1A是正定的,但通常不是对称的,因此我们不能直接对(5.6)实施CG方法.此时我们需要定义一个新的内积:P-内积,即(x,y)P≜(P x,y)=(x,P y).(5.7)由于P是对称正定的,所以这个内积定义是有效的.在该内积下,有(P−1Ax,y)P=(Ax,y)=(x,Ay)=(x,P(P−1Ay))=(x,P−1Ay)P,即P−1A关于P-内积是自伴随的,即在P-内积意义下,P−1A是“对称”的.这时可以用CG 方法来求解预处理方程组(5.6),但需要将欧拉内积改为P -内积.引入辅助变量z k ≜P −1r k ,则求解预处理方程组(5.6)的CG 方法可描述如下:算法1.3基于P -内积的CG 方法1:给定初值x (0),(相对)精度要求ε>0和最大迭代步数IterMax 2:计算r 0=b −Ax (0)和β=∥r 0∥23:令z 0=P −1r 0,p 1=z 04:for k =1to IterMax do 5:ξk =(z k −1,z k −1)P (P −1Ap k ,p k )P =(r k −1,z k −1)(Ap k ,p k )6:x (k )=x (k −1)+ξk p k 7:r k =r k −1−ξk Ap k 8:relres =∥r k ∥2/β9:if relres <εthen10:break 11:end if 12:z k =z k −1−αk P −1Ap k =P −1r k 13:µk =(z k ,z k )P (z k −1,z k −1)P =(r k ,z k )(r k −1,z k −1)14:p k +1=z k +µk p k 15:end for1.3预处理GMRES方法左预处理GMRES对于非对称情形的GMRES方法,也有三种预处理方式.但与对称正定情形不同的是,这三种方式并不等价,而且有时效果会相差很大.首先考虑左预处理方式,即P−1Ax=P−1b(5.8)记r k=b−Ax(k)为原问题的残量,则预处理方程组(5.8)的残量为˜r k=P−1r k将GMRES用于求解方程组(5.8),即可得左预处理GMRES方法.算法1.4左预处理GMRES方法1:给定初值x(0)和(相对)精度要求ε>02:计算˜r0=P−1(b−Ax(0))和β=∥˜r0∥23:令v1=˜r0/β,ξ=βe14:for j=1,2,...do5:˜w j=Av j6:w j=P−1˜w j%Apply preconditioning,solve P w j=˜w j for w j7:for i=1,2,...,j do8:h ij=(w j,v i)9:w j=w j−h ij v i10:end for11:h j+1,j=∥w j∥212:for i=1,2,...,j−1do%Apply G j−1,...,G1to last column of H j+1,j13:[h i,jh i +1,j ]=[c i s i −s i c i ][h i,j h i +1,j ]14:end for 15:if h j +1,j =0then 16:set m =j and break 17:end if 18:v j +1=w j /h j +1,j 19:if |h j,j |>|h j +1,j |then %Form the Givens rotation G j 20:c j =1√1+τ2,s j =c j τwhere τ=h j +1,j h j,j 21:else 22:s j =1√1+τ2,c j =s j τwhere τ=h j,j h j +1,j 23:end if 24:h j,j =c j h j,j +s j h j +1,j ,h j +1,j =0%Apply G j to last column of H j +1,j25:[ξjξj+1]=[c j s j−s j c j][ξj]%Apply G j to the right-hand side26:if|ξj+1|/β<εthen%Check convergence:∥˜r j∥2=|ξj+1|27:set m=j and break28:end if29:end for30:计算y(m)=R−1mξ(m),其中R m=H(1:m,1:m),ξ(m)=ξ(1:m)31:计算近似解x(m)=x(0)+V m y(m)需要指出的是,在左预处理GMRES方法中,停机准则中的|ξj+1|并不是原方程组残量r k的模,而是预处理残量˜r k=P−1(b−Ax(k))的模.如果要计算真实残量,除了显式计算外,没有其它更好的办法.所以,如果停机准则是基于真实残量的,则需要额外去计算r k,这就会增加很多额外的运算量.右预处理GMRES如果用P做右预处理,则对应的预处理方程组为AP−1u=b,x=P−1u.(5.9)给定迭代初值x(0),则u(0)=P x(0).相应的预处理初始残量为˜r0=b−AP−1u(0)=b−Ax(0)=r(0),即预处理残量与原方程组残量是一样的.因此无需计算u(0).设u(m)=u(0)+V m y m 是迭代m步后的近似解,则对应的原方程组近似解x(m)为x(m)=P−1u(m)=P−1(u(0)+V m y m)=x(0)+P−1V m y m,原方程组的残量为r m=b−Ax(m)=r0−AP−1V m y m=βv1−V m+1H m+1,m y m=V m+1(βe1−H m+1,m y m).由于y m是最小二乘问题∥βe1−H m+1,m y∥2miny∈R m的解,因此∥r m∥2=∥βe1−H m+1,m y m∥2=β·|q1(m+1)|.这与不带预处理的GMRES方法是一样的.因此在实际求解过程中我们可以直接得到原方程组残量的模,而无需计算u(m)和x(m).这是右预处理方式与左预处理方式的主要区别.右预处理GMRES方法具体描述如下:算法1.5右预处理GMRES方法1:给定初值x(0)和(相对)精度要求ε>0 2:计算r0=b−Ax(0)和β=∥r0∥23:令v1=˜r0/β,ξ=βe14:for j=1,2,...do5:˜w j=P−1v j%Apply preconditioning 6:w j=A˜w j7:for i=1,2,...,j do8:h ij=(w j,v i)9:w j=w j−h ij v i10:end for11:h j+1,j=∥w j∥212:for i =1,2,...,j −1do %Apply G j −1,...,G 1to the last column of H j +1,j 13:[h i,j h i +1,j ]=[c i s i −s i c i ][h i,jh i +1,j]14:end for15:if h j +1,j =0then16:set m =j and break17:end if18:v j +1=w j /h j +1,j19:if |h j,j |>|h j +1,j |then %Form the Givens rotation G j 20:c j =1√1+τ2,s j =c j τwhere τ=h j +1,jh j,j21:else22:s j =1√1+τ2,c j =s j τwhere τ=h j,jh j +1,j23:end if24:h j,j =c j h j,j +s j h j +1,j ,h j +1,j =0%Apply G j to last column of H j +1,j 25:[ξj ξj +1]=[c j s j −s j c j ][ξj 0]%Apply G j to the right-hand side 26:if |ξj +1|/β<εthen %Check convergence:∥r j ∥2=|ξj +1|27:set m =j and break 28:end if29:end for30:计算y (m )=R −1m ξ(m ),其中R m =H (1:m,1:m ),ξ(m )=ξ(1:m )31:计算近似解x (m )=x (0)+P −1V m y (m )两边预处理方式如果预处理子P是由乘积形式给出的,即P=P L P R,则可以构造下面的两边预处理方程组P−1 L AP−1Ru=P−1Lb,x=P−1Ru.将GMRES用于求解上述线性方程组,则可得到两边预处理GMRES方法.具体推导过程与前面类似,留作练习.需要指出的是,对于两边预处理GMRES方法,其预处理后的残量为˜r k=P−1L(b−Ax(k)).1.4左、右预处理方式的最优性质设x(0)是迭代初始值,记x(k)L 和x(k)R分别是左预处理GMRES方法和右预处理GM-RES方法迭代k步后得到的近似解.根据GMRES方法的最优性质可知,对于左预处理GMRES方法,有x(k)L=argminx∈x(0)+K k(P−1A,P−1r0)∥P−1(b−Ax)∥2(5.10)而对于右预处理GMRES,我们有x(k)R=P−1u(k),其中u(k)=argminu∈u(0)+K k(AP−1,r0)∥b−AP−1u∥2.通过变量代换x=P−1u,可得x(k)R=argminx∈x(0)+P−1K k(AP−1,r0)∥b−Ax∥2(5.11)又P−1K k(AP−1,r0)=P−1span{r0,AP−1r0,(AP−1)2r0,...,(AP−1)k−1r0}=span{P−1r0,P−1AP−1r0,P−1(AP−1)2r0,...,P−1(AP−1)k−1r0}=span{P−1r0,(P−1A)P−1r0,(P−1A)2P−1r0,...,(P−1A)k−1P−1r0} =K k(P−1A,P−1r0),由(5.10)和(5.11)可知,x(k)L 和x(k)R属于同一个子空间中.不同之处在于,x(k)L极小化的是预处理后残量的范数,即∥P−1(b−Ax)∥2,而x(k)R极小化的是原方程组残量的范数,即∥b−Ax∥2.1.5Flexible GMRES方法在前面的讨论中,我们假定预处理子P是固定不变的.事实上,我们可以在每步迭代中取不同的预处理子.这就是Flexible GMRES(FGMRES)方法的基本思想.下面介绍右预处理方式的Flexible GMRES方法.我们知道,在右预处理GMRES方法1.5中,近似解可以写成下面的形式:x(m)=x(0)+P−1V m˜y=x(0)+m∑j=1(P−1v j)˜y j.(5.12)这意味着x(m)−x(0)是向量P−1v1,P−1v2,...,P−1v m的线性组合.而这些向量可由算法1.5中的第5步计算得到.但事实上,我们无需存储这些向量,因为它们是由同一个预处理子P作用到向量v′j s上所得到的.因此,在计算近似解x(m)时,只需将P作用到v1,v2,...,v m的线性组合上即可,即V m˜y上.但如果每次迭代时采用不同的预处理子,即将算法1.5中的第5步改为5':solve P j w=v j for w,that is,w=P−1j v j因此,我们很自然地想到用下面的方式来获取近似解:x(m)=x(0)+m∑j=1(P−1jv j)˜y j=x(0)+Z m˜y,其中Z m=[z1,z2,...,z m].这z j=P−1jv j,向量˜y则跟之前一样计算.需要注意的是,此时我们需要存储所有的向量z j=P−1jv j,j=1,2,...,m.这就是FGMRES与右预处理GMRES方法的根本区别.FGMRES算法描述如下:算法1.6Flexible GMRES(FGMRES)1:给定初值x(0)和(相对)精度要求ε>02:计算r0=b−Ax(0)和β=∥r0∥23:setξ=βe1and v1=r0/β4:for j=1,2,...do5:solve P j z j=v j for z j,that is,z j=P−1v j%Apply the preconditioner 6:w j=Az j7:for i=1,2,...,j do%Arnoldi process8:h ij=(w j,v i)9:w j=w j−h ij v i10:end for11:h j+1,j=∥w j∥212:if h j +1,j =0then 13:set m =j and break 14:end if15:v j +1=w j /h j +1,j16:for i =1,2,...,j −1do %Apply G j −1,...,G 1to the last column of H i +1,j17:[h i,j h i +1,j ]=[c i s i −s i c i ][h i,jh i +1,j ]18:end for19:if |h j,j |>|h j +1,j |then %Form the Givens rotation G j20:set c j =1√1+τ2,s j =c j τwhere τ=h j +1,jh j,j 21:else 22:set s j =1√1+τ2,c j =s j τwhere τ=h j,j h j +1,j23:end if24:h j,j =c j h j,j +s j h j +1,j ,h j +1,j =0%Apply G j to last column of H j +1,j 25:[ξj ξj +1]=[c j s j −s j c j ][ξj0]%Apply G j to the right-hand side 26:if ∥ξj +1∥2<εthen %Check convergence:∥˜r ∥2=ξj +127:set m =j and break28:end if 29:end for30:compute ˜y =R −1m ξm with R m =H (1:m,1:m )upper-triangular and ξm =ξ(1:m )31:compute the approximate solution x (m )=x (0)+Z m ˜y1.6预处理BiCGSTAB方法右预处理BiCGSTAB与带预处理的GMRES方法类似,为了使得计算过程中得到的残量是真实残量(而不是预处理后的残量),我们一般采用右预处理方式,即AM−1u=b with u=Mx,(5.13)其中M为预处理子.在BiCGSTAB方法中,主要迭代过程为:q k=r k−1−αk Ap k withαk=(r k−1,˜r0) (Ap k,˜r0),r k=q k−ωk Aq k withωk=(q k,Aq k) (Aq k,Aq k),x(k)=x(k−1)+αk p k+ωk q k,p k+1=r k+βk(p k−ωk Ap k)withβk=αkωk·(r k,˜r0)(r k−1,˜r0),其中k=1,2,...,p1=r0.记˜p k=M−1p k,˜q k=M−1q k,则针对预处理后的线性方程组(5.13)的BiCGSTAB迭代过程为:q k=r k−1−αk A˜p k withαk=(r k−1,˜r0) (A˜p k,˜r0),r k=q k−ωk A˜q k withωk=(q k,A˜q k) (A˜q k,A˜q k),x(k)=x(k−1)+αk˜p k+ωk˜q k,p k+1=r k+βk(p k−ωk A˜p k)withβk=αkωk·(r k,˜r0)(r k−1,˜r0).整理后即可得下面的预处理BiCGSTAB方法:算法1.7Preconditioned BiCGSTAB1:Given an initial guess x(0)and stopping toleranceε>02:compute r0=b−Ax(0)3:choose˜r0such that(r0,˜r0)=0%one common choice is˜r0=r0 4:set p1=r05:for k=1,2,...do6:ρk−1=(r k−1,˜r0)%ifρk−1=0,then breakdown7:if k>1then8:βk−1=(αk−1/ωk−1)·(ρk−1/ρk−2)9:p k=r k−1+βk−1(p k−1−ωk−1v k−1)10:end if11:˜p k=M−1p k%solve M˜p k=p k for˜p k12:v k=A˜p k13:αk=ρk−1/(v k,˜r0)%if(v k,˜r0)=0,then breakdown 14:q k=r k−1−αk v k15:x(k)=x(k−1)+αk˜p k%if∥q k∥<ε,then converged and stop 16:˜q k=M−1q k%solve M˜q k=q k for˜q k17:w k=A˜q k18:ωk=(q k,w k)/(w k,w k)%ifωk=0,then breakdown19:r k=q k−ωk w k20:x(k)=x(k)+ωk˜q k21:if∥r k∥2<εthen%check convergence22:stop%converged23:end if24:end for2基于矩阵分裂的预处理技术Finding a good preconditioner to solve a given sparse linear system is oftenviewed as a combination of art and science.Theoretical results are rare andsome methods work surprisingly well,often despite expectations.—Saad,2003.预处理能否取得成功的关键是能否找到一个好的预处理子.预处理子的构造与问题本身是密切相关的,通用的“好预处理子”是不存在的.关于预处理技术的理论分析很少,大多数情况下只能根据经验来构造.尽管如此,在实际应用中,这些根据经验构造出来的预处理子往往能取得很好的效果,有时效果之好会大大出乎人们的意料之外.一般来说,预处理子可以分为两大类:(a)代数预处理子(Algebraic Preconditioner),即仅仅根据所给方程组的系数矩阵来构造预处理子.(b)专用预处理子(Problem-Specific Preconditioner),即根据问题的物理背景所构造的专用预处理子.显然,由于专用预处理子充分利用了问题的物理背景知识,所以它们往往具有很好的数值表现,如多重网格,区域分解,快速变换等等.但它们严重依赖于原问题的物理背景,因此不具有通用性.我们这里只介绍代数预处理子,即仅仅根据所给的系数矩阵来构造预处理方法.这种预处理方法具有一定的通用性.考虑线性方程组Ax=b,A∈R n×n,对A做如下的矩阵分裂:A=M−N(5.14)其中M非奇异,则可以得到下面的迭代方法x(k+1)=M−1Nx(k)+M−1b=x(k)+(M−1b−M−1A)x(k).k=0,1,2,....这等价于求解下面的方程组M−1Ax=M−1b.(5.15)这就是与矩阵分裂(5.14)相对应的预处理线性方程组.。
第八讲预处理技术
5/20
预处理方式的选择
若 A 对称正定, 用 CG 求解, 这三种方式的预处理效果基本一致. 但若 A 非对称 (特别是非正规) 情形, 效果可能会相差很大. 实际使用中, 选取哪种预处理方式, 根据问题本身和所用的算法来确定. 对于 GMRES, 右预处理比较合适, 因为预处理后的残量就是真实残量. † 需要指出的是, 在实际求解预处理后的方程组时, 我们并不需要显式 地计算 P −1 (除非 P −1 非常容易计算), 更不需要显式地计算 P −1 A.
1: 5: 6: 7: 8: 9: 10:
x(k) = x(k−1) + ξk pk rk = rk−1 − ξk Apk zk = zk−1 − αk P −1 Apk = P −1 rk (rk , zk ) ( zk , z k ) P µk = = ( z k −1 , z k −1 ) P (rk−1 , zk−1 ) pk+1 = zk + µk pk end for
14/20
右预处理
AP −1 u = b, x = P −1 u.
给定迭代初值 x0 , 则 u(0) = P x(0) . 相应的预处理初始残量为 r ˜0 = b − AP −1 u(0) = b − Ax(0) = r(0) . 因此无需计算 u(0) . 设 u(m) = u(0) + Vm ym 是迭代 m 步后的近似解, 则 x(m) = P −1 u(m) = P −1 (u(0) + Vm ym ) = x(0) + P −1 Vm ym , 真实残量为 rm = b − Ax(m) = r0 − AP −1 Vm ym = Vm+1 (βe1 − Hm+1,m ym ). 因此 ∥rm ∥2 = ∥βe1 − Hm+1,m ym ∥2 = β · |q1 (m + 1)|.
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
预处理计算方法
●多介质过滤器、活性炭过滤器、保安过滤器和树脂罐型号的选择:
主要看RO终端出水流量,以RO设备对水的利用率计算预处理的流量.一般小型RO对水利用率大概为50---60%,大型系统在75%，最高不超过80%。
举例说明相关问题：
假设一级RO产水30 m3/h,回收率取75%，预处理产水量40 m3/h，暂不考虑预处理预备自用水率；
◆多介质过滤器产水量40 m3/h，设计流速取8---10m/h，根据Q=SV
【Q流量、S罐体截面积、V设计流速,S=π R2.罐体直径d=2R】
可以计算出罐体的直径，高度查厂家手册。
◆活性炭过滤器产水量40 m3/h，设计流速取10---12m/h，根据Q=SV
【Q流量、S罐体截面积、V设计流速,S=π R2.罐体直径d=2R】可
以计算出罐体的直径，高度查厂家手册。
◆软化器产水量40 m3/h，设计流速取15---20m/h，根据Q=SV【Q流
量、S罐体截面积、V设计流速,S=π R2.罐体直径d=2R】可以计算
出罐体的直径，高度查厂家手册。
◆保安过滤器产水量40 m3/h，滤芯通量取10 m3/ m3•h，按照滤芯的
规格，40〞滤芯外表面积是0.2 m2，即单只滤芯产水2 m3/h,（滤
芯的外径基本是63㎜，长度取1m，计算下来的表面积是0.1979 m2）
按照这个通量计算下来就是2m3/40〞，滤芯总数20支，选择20芯
40〞的过滤器即可。