CN109492256B - 一种动态系统状态概率映射矩阵多向搜索方法 - Google Patents

Abstract

本发明公开了一种动态系统状态概率映射矩阵多向搜索方法，包括以下内容：建立系统仿真和故障注入模型，系统状态空间离散化，系统始发栅元点随机抽样，配置系统部件状态，设置仿真时间步长、仿真测试次数，通过函数赋值将系统始发栅元点信息导入仿真模型中，执行并完成仿真测试，统计系统轨迹在栅元中的迁移分布，计算系统状态条件迁移概率，生成系统状态转移概率矩阵，建立系统动态概率映射网络；进行搜索参数初始化，执行矩阵动态搜索分析，实现系统重要风险序列路径识别和查找，计算并预测系统可靠度随时间的概率分布。本发明可用于系统动态轨迹及其状态转移概率矩阵生成，支持矩阵多向搜索计算，获得动态可靠性评价指标。

Description

一种动态系统状态概率映射矩阵多向搜索方法

技术领域

本发明涉及数理统计和仿真分析领域，具体涉及一种动态系统状态概率映射矩阵多向搜索方法。

背景技术

近年来，随着数字化系统技术的高速发展，核电厂也正面临其仪表与控制系统数字化升级换代及全面应用的趋势。仪表与控制系统作为核电厂的大脑神经，牵带着系统中每一参数的监督和控制，一旦发生系统故障或失效，将直接影响或损害核电机组的安全运行。而由于核电具有的潜在放射性物质释放风险，核电厂对数字化技术产品的安全可靠性标准也有区别于其他行业的特殊的、严格的要求。

数字化技术产品相比于传统的模拟式仪表和控制器件应用优势明显，但其独特的技术特点(如软件和硬件的集成，交互耦合失效等)也引入了新的挑战，加大了核电厂数字化仪表和控制系统的认证难度。传统的可靠性建模分析方法如故障树/事件树方法已不能满足数字化仪表和控制系统复杂的动态交互过程分析需求，国内外相关领域专家也都在致力于寻求一种动态的可包罗数字化仪表与控制系统特有复杂特性的新型可靠性分析方法，进而实现对数字化过程交互的动态可靠性和安全评价，并为核电厂数字化仪表和控制系统授权和应用进程提供审查技术支持。

发明内容

本发明的目的是针对现有技术的不足，提出了一种动态系统状态概率映射矩阵多向搜索方法，所述方法提供了一套可以通过故障注入的系统仿真模型，模拟系统状态在不同离散化状态空间栅元中的动态迁移过程，以此估算并建立系统状态概率映射关系及转移概率矩阵，并通过搜索算法实现系统关键序列路径生成和可靠性定量评估。所述方法能够应用于核电厂数字化仪表与过程控制系统中，针对其动态交互耦合特性，实现过程控制系统故障注入模拟测试、以及自动实现系统可靠性建模更新和评价分析的动态概率映射矩阵搜索，实现对数字化过程交互的动态可靠性和安全评价。

本发明的目的可以通过如下技术方案实现：

一种动态系统状态概率映射矩阵多向搜索方法，所述方法包括以下内容：

对数字化系统部件进行失效模式、失效机理及影响分析，定义系统控制单元物理状态，建立表述过程控制系统部件状态转移的有限状态机模型；

依据过程控制系统控制法则建立其仿真模型，并根据数字化系统部件失效机理建立故障注入模型；

定义过程控制系统状态空间，对系统连续过程变量进行空间栅元划分，并通过随机函数抽样栅元中系统状态的始发点；

设置系统状态配置参数、仿真测试次数、仿真时长，并通过函数赋值在仿真模型中导入随机抽样系统状态点的位置；

对系统进行仿真测试，针对每一系统始发栅元点，统计其系统运行轨迹落入不同栅元中的点的个数，计算系统状态条件转移概率，生成系统状态转移概率矩阵，建立系统动态概率演变模型；

设置系统初始状态

搜索深度n、系统初始状态概率P和截断准则ε，其中Ω表示过程控制系统状态空间，V_j表示将过程控制系统状态空间Ω划分的互不重叠的离散化栅元；进行多向矩阵动态搜索分析，通过前向搜索实现系统状态演变重要序列路径生成，以及通过回溯搜索实现系统失效重要风险路径识别；最后计算和预测系统状态在不同时刻的可靠度。

进一步地，所述定义过程控制系统状态空间，对系统连续过程变量进行空间栅元划分，并通过随机函数抽样栅元中系统状态的始发点过程中，系统状态时空离散化方法为：

将由L个连续变量构成的L维过程控制系统状态空间Ω划分为互不重叠的离散化栅元V_j＝{x_l:a_l,j≤x_l≤b_l,j；j＝1,...,J_l；l＝1,...,L}，其中x_l表示具体的某一空间栅元，a_l,j和b_l,j分别表示该空间栅元的下边界值和上边界值，J_l表示针对每一变量维度所划分的区间数；系统动态能够看作是特定系统物理控制部件状态组合m(m＝1,...,M)下系统状态在离散化空间栅元V_j中随时间的概率转移过程，系统状态转移时间步长使用Δt表示，空间栅元划分遵循最优原则，应选择适当的栅元尺寸和时间步长以确保系统状态转移尽量发生在相邻栅元之间；假设部件状态配置，即数字化部件状态的组合m变化只可能发生在t+Δt时刻，而不会出现在[t,t+Δt)时间段内，则系统在t+Δt和t时刻的状态前向相关性能够通过Pⁿ⁺¹＝Q·Pⁿ表示；如果是反向搜索应用，则系统状态之间的相关性能够通过Pⁿ＝[Q^TQ]^{- 1}Q^T·Pⁿ⁺¹表示；其中，列向量Pⁿ＝[pⁿ(1,1),pⁿ(1,2),...,pⁿ(1,J),...,pⁿ(M,1),...,pⁿ(M,J)]^T表示系统状态在特定搜索深度n下的概率；系统状态转移概率Q＝q(m,j|m',j',Δt)；m,m'＝1,...,M；j,j'＝1,...,J由q(m,j|m',j',Δt)＝g(j|j',m',Δt)·h(m|m',j'→j,Δt)计算得到；其中，h(m|m',j'→j,Δt)表示部件状态配置在系统过程变量栅元j'→j转移过程中由m'变为m的概率；h(m|m',j'→j,Δt)为系统控制单元随机状态变化的全概率映射；假设部件失效之间相互独立，h(m|m',j'→j,Δt)能够简单地表示成每一单独控制单元失效概率的乘积；g(j|j',m',Δt)表示在给定部件状态配置m'下和时间步长Δt内，系统过程变量由栅元j'进入栅元j的条件概率；g(j|j',m',Δt)矩阵数据值能够通过系统控制数学模型或仿真模型/模拟器获得；通过对系统状态转移概率矩阵Q的计算，实现对系统状态概率映射关系模型的构建；

所述系统状态的始发点随机抽样方法为：通过一随机抽样函数每次从指定栅元V_j中随机抽取一个或多个数值作为系统状态的位置起始点用于系统轨迹模拟生成，随机抽样函数服从均匀分布；系统测试时，以栅元中的点代表栅元。

进一步地，所述g(j|j',m',Δt)矩阵数据值能够通过系统控制数学模型采用连续积分法或仿真模型/模拟器采用等权点积分获得；此处的g(j|j',m',Δt)元素值通过仿真模型模拟得到。

进一步地，所述系统仿真测试通过命令程序语句实现，测试结果包括系统落入栅元点的位置、数字化部件故障模式的符号参数以及最大和最小失效指数，并将所述测试结果存入数据库中。

进一步地，所述系统状态条件转移概率通过以下等权点积分方法近似得到：假若系统测试过程中在栅元j'中随机抽取了u_j'个数作为系统状态起始位置点，这u_j'个数均匀分布在栅元j'中，经系统状态发展轨迹模拟，如果在时间步长Δt内有w_j个系统轨迹点落入栅元j，则系统状态条件转移概率矩阵近似于

进一步地，所述多向矩阵动态搜索分析以系统指定初始状态s为起点，通过搜索系统状态转移概率矩阵Q中的非零元素，并将观测到的非零元素作为父节点继续下一步搜索，迭代循环直至指定搜索深度n；前、后搜索同步进行，每一搜索步都涉及符合元素的状态转移概率矩阵Q的更新计算，实行广度搜索优先，以遍历所有可能的序列路径；

其中，系统状态演变序列路径通过不同搜索步上系统状态栅元点前后相连构成的事件序列线路生成，系统状态演变序列路径的概率为每一搜索步分支概率的乘积；系统重要序列路径筛选方法如下：系统状态序列路径的重要性体现在其概率值排序上，通过截断准则ε以筛选出系统状态发展重要序列路径以及系统失效重要风险路径。

进一步地，所述系统状态在不同时刻的可靠度计算方法如下：特定搜索深度n下，对系统状态进行分组归类，属于同一系统状态的序列分支归为一组，系统状态在n·Δt时刻的可靠度为组内所有序列分支概率之和，系统状态概率的预测通过Pⁿ⁺¹＝Q·Pⁿ实现。

本发明与现有技术相比，具有如下优点和有益效果：

1、本发明提供的动态系统状态概率映射矩阵多向搜索方法，由系统仿真模型及故障注入测试生成的系统轨迹及系统状态转移概率矩阵，考虑了系统随机失效、时序关系、非线性特性、偶然不确定性和知不确定性，更符合系统动态过程，仿真意义更加贴合实际。

2、本发明提供的动态系统状态概率映射矩阵多向搜索方法，可以通过增加系统仿真测试的次数，自动更新系统状态转移概率映射矩阵模型；且系统状态转移概率矩阵可支持多方向搜索，并从概率风险的角度筛选出系统重要序列路径，识别系统设计和运行过程中的薄弱环节。

附图说明

图1为本发明实施例中数字化部件故障注入通选模型。

图2为本发明实施例中基于模型测试的系统状态转移概率矩阵计算程序框图。

图3为本发明实施例中系统状态转移概率矩阵通用搜索计算程序框图。

图4为本发明实施例中系统运行动态概率预测曲线图。

图5为本发明实施例中系统轨迹发展示意图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例：

本实施例提供了一种动态系统状态概率映射矩阵多向搜索方法，属于数理统计和仿真分析领域，主要是模拟复杂过程控制系统轨迹生成，实现动态系统可靠性建模和矩阵搜索分析。本实施例所述方法的算法程序和软件应用平台分别是以MATLAB/SIMULINK 2016和Eclipse 4.6平台编写的，其主要功能为：辅以故障模式及影响分析，对数字化系统部件故障模式和失效机理进行分类，建立有限状态机模型，模拟数字化系统部件随机状态转移；根据数字化系统部件的失效机理构建系统仿真和故障注入模型，对系统连续过程变量状态空间进行离散化；通过栅元点随机抽样和系统轨迹生成模拟测试，建立系统状态转移概率映射矩阵模型，设定系统模型分析用参数输入，实现特定搜索深度下系统重要风险路径查找和识别，定量预测系统可靠性动态变化，最终完成对大型复杂过程系统动态可靠性和安全评价和验证。

整个软件应用平台包括可视化图形显示界面、图形化建模单元、系统状态空间栅元定义制表、数据导入导出功能、系统参数初始化功能、系统仿真和故障注入测试功能、系统状态转移概率映射矩阵模型生成和更新、通用矩阵搜索算法、系统动态可靠性分析功能等几个模块。

进入系统应用平台后，导入系统数字化部件故障模式及影响分析结果，绘制系统部件状态转移有限状态机模型，计算部件状态配置概率转移矩阵h(m|m',j'→j,Δt)。与此同时，对系统状态时间和空间进行离散化定义，搭建维度化栅元结构框架，由系统仿真和故障注入测试功能单元设置系统起始栅元点参数、部件状态配置及测试试验次数，然后通过功能按钮弹出“模拟计算开始”对话框，确认进入系统动态测试环节，计算和动态更新系统状态条件转移概率g(j|j',m',Δt)，待模拟达到指定次数要求，弹出“模拟计算结束”对话框，计算结果显示于对话框中，结束本次仿真测试调用。重新设置新的系统起始栅元点参数、部件状态配置等信息，重复上述模拟计算过程，直至系统状态条件转移概率矩阵覆盖所有可能系统栅元点取样。仿真测试完成之后，利用矩阵乘积运算得到系统状态转移概率映射关系矩阵Q。

本实施例可以基于系统状态转移概率矩阵Q进行多向同步搜索推理分析。搜索流程为：首先通过软件应用平台进行算法初始化参数设置，包括指定系统初始状态s及其概率P、搜索深度n和截断准则ε，待初始化参数设置完成后，确定并开始矩阵搜索，搜索过程中，被观测到的目标矩阵元素将用作父节点开始新的子节点搜索，期间与目标元素相关联的转移分支将对矩阵Q更新计算。如图3所示，算法会对每一搜索步上的所有可能序列路径进行概率计算，并通过截断准则对低概率事件序列分支进行剔除，筛选出系统失效重要风险路径并存入数据库中，然后对保留下来的系统重要序列分支进行概率更新，向前和向后同时更新系统循环变量n，进入下一次迭代循环。矩阵搜索过程中，被观测到的重要序列分支数目会伴随搜索步深度在搜索结果对话框中滚动更新，每一步得到序列路径总数会在搜索步结束之后存入数据库中以进行统计分析。

系统状态概率定量计算也是针对具体搜索步进行，首先对特定系统状态所在序列分支进行归组，然后对组内所有进入目标系统状态的序列分支概率求和，累加数值即为特定系统状态当下的概率值。

系统搜索分析计算完成后，将系统失效重要风险路径以向量的形式存入数据库表格中，用于数据处理分析和备份，基于分析结果数据绘制系统状态概率时间关系曲线。

接下来，通过对核电厂数字化仪表与过程控制系统动态建模仿真和可靠性分析，验证本实施例所述方法的有效性，包括以下步骤：

第一步：针对示例系统中数字化控制单元进行故障模式、故障机理及影响分析，定义系统控制单元物理状态。系统由两个进水控制器、一个出水控制器组成，每一控制单元包含运行、失效“关”和失效“开”三种状态，状态转移发生在运行和失效之间。

第二步，建立水位控制示例系统MATLAB/SIMULINK仿真模型，并针对数字化部件建立其故障注入模型，如图1所示。

第三步：对水位过程变量进行区间栅元划分，系统水位有效控制区域为-3米≤x≤3米，控制区域内栅元长度为2米，定义后的水位状态栅元分别为：x₁<-3米，-3米≤x₂≤-1米，-1米<x₃≤1米，1米<x₄≤3米，x₅>3米。

第四步：结合第一步给出系统部件状态定义和第三步给出的系统过程变量状态划分，得出系统控制部件状态组合数为M＝3³＝27，离散化后的系统过程变量状态空间栅元数为J＝5，计算得出系统状态总数为N(s)＝M·J＝3×3×3×5＝135。同时，取仿真时间步长Δt＝30分钟。

第五步：配置系统部件组合状态m，将数字化部件失效影响程度和始发栅元点位置分别使用两个不同的随机抽样函数表示，抽样时间步长与仿真时间步长一致。本实施例中，系统部件状态数字编码为一3维行向量，控制部件三种状态模式定义分别使用1、2、3数字编码表示，并将所有的27组部件配置向量存入元胞数组中。使用循环变量m(m＝1,2,...,27)逐个取出元胞数组中的数字编码进行部件状态配置。

数字化部件失效影响程度使用随机信号生成器Amount_{In/Decrease％}＝normrnd(μ,σ,1,u_j)描述，生成的随机信号服从正态分布，其中μ＝0，σ＝0.5，u_j＝100。控制部件失效影响程度表示部件失效引起的相对于其正常输出而产生的最大增量和减量百分比，最大增量和减量分别对应于正态分布函数的右边正值取数和左边负值取数，最大故障输出相对百分比为1，通过限值函数将正态分布取数中绝对值超出1的数值部分截断，保留最大±1数值。实施过程中，通过绝对值正负跳转形式实现部件失效影响程度和部件运行/失效模式的匹配和统一，如部件失效开对应正输出，则所有的负抽样数值自动取其正绝对值作为附加失效输出。

随后，编写随机数生成函数p_j-departure＝a_j+(b_j-a_j).*rand(u_j)用于抽样系统始发栅元点。其中，p_j-departure表示针对栅元j选取的始发点位置，b_j和a_j分别表示栅元j的上、下边界值，rand(u_j)表示从开区间(0,1)随机抽取u_j个系统始发栅元点，栅元抽样点服从均匀分布。本实施例中，对27组部件配置中的每一种情况执行100次仿真测试，测试数据存入数据库中，见图2。

第六步：根据系统仿真输出结果判断系统根轨迹终点位置所处栅元，统计每一栅元中落入点的个数w_j，将栅元中落入点的数目w_j除以始发栅元点数目u_j，得到对应栅元条件迁移概率，并进一步通过等权点积分计算式求得所有系统状态条件转移概率矩阵g(j|j',m',Δt)。

第七步：由系统部件失效率和修复率及仿真时间步长计算部件概率转移矩阵h(m|m',j'→j,Δt)，通过矩阵乘积运算得到系统状态转移概率矩阵Q。

第八步：对搜索算法进行参数初始化，设定系统起始状态栅元为s∈V₃，搜索深度n＝8、系统初始状态概率P＝1、截断准则ε＝0。

第九步：依据第八步给出的初始化参数输入，执行并完成多向搜索分析，将每一层次搜索结果存入数据库中，算法实施过程见图3。

第十步：针对特定的目标系统状态s∈V₁，归组所有导致系统低水位失效的序列分支，并进行概率求和，得到系统低水位失效概率随时间函数分布，结合图表在应用平台界面上显示，如图4所示。同样，在系统初始状态和终点状态一致的前提下，对部件配置变化次序相同的序列路径进行归组，通过求和计算出系统动态发展情景的概率分布。依据序列或情景分支概率权重值大小，从风险的角度对所有路径分支进行重要性排序，然后基于系统失效重要序列情景，绘制系统轨迹图，如图5所示。

以上所述，仅为本发明专利较佳的实施例，但本发明专利的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明专利所公开的范围内，根据本发明专利的技术方案及其发明专利构思加以等同替换或改变，都属于本发明专利的保护范围。

Claims (7)

1.一种动态系统状态概率映射矩阵多向搜索方法，其特征在于，所述方法包括以下内容：

对数字化系统部件进行失效模式、失效机理及影响分析，定义系统控制单元物理状态，建立表述过程控制系统部件状态转移的有限状态机模型；

依据过程控制系统控制法则建立其仿真模型，并根据数字化系统部件失效机理建立故障注入模型；

定义过程控制系统状态空间，对系统连续过程变量进行状态空间栅元划分，并通过随机函数抽样栅元中系统状态的始发点；

设置系统状态配置参数、仿真测试次数、仿真时长，并通过函数赋值在仿真模型中导入随机抽样系统状态始发点的位置；

对系统进行仿真测试，针对每一系统始发栅元点，统计其系统运行轨迹落入不同栅元中的点的个数，计算系统状态条件转移概率，生成系统状态转移概率矩阵，建立系统动态概率演变模型；

设置系统初始状态

搜索深度n、系统初始状态概率P和截断准则ε，其中Ω表示过程控制系统状态空间，V_j表示将过程控制系统状态空间Ω划分的互不重叠的离散化栅元；进行多向矩阵动态搜索分析，通过前向搜索实现系统状态演变重要序列路径生成，以及通过回溯搜索实现系统失效重要风险路径识别；最后计算和预测系统状态在不同时刻的可靠度。

2.根据权利要求1所述的一种动态系统状态概率映射矩阵多向搜索方法，其特征在于，所述定义过程控制系统状态空间，对系统连续过程变量进行空间栅元划分，并通过随机函数抽样栅元中系统状态的始发点过程中，系统状态时空离散化方法为：

将由L个连续变量构成的L维过程控制系统状态空间Ω划分为互不重叠的离散化栅元V_j＝{x_l:a_l,j≤x_l≤b_l,j；j＝1,...,J_l；l＝1,...,L}，其中x_l表示具体的某一空间栅元，a_l,j和b_l,j分别表示该空间栅元的下边界值和上边界值，J_l表示针对每一变量维度所划分的区间数；系统动态能够看作是特定系统物理控制部件状态组合m(m＝1,...,M)下系统状态在离散化空间栅元V_j中随时间的概率转移过程，系统状态转移时间步长使用Δt表示，空间栅元划分遵循最优原则，应选择适当的栅元尺寸和时间步长以确保系统状态转移尽量发生在相邻栅元之间；假设部件状态配置，即数字化部件状态的组合m变化只可能发生在t+Δt时刻，而不会出现在[t,t+Δt)时间段内，则系统在t+Δt和t时刻的状态前向相关性能够通过Pⁿ⁺¹＝Q·Pⁿ表示；如果是反向搜索应用，则系统状态之间的相关性能够通过Pⁿ＝[Q^TQ]^-1Q^T·Pⁿ⁺¹表示；其中，列向量Pⁿ＝[pⁿ(1,1),pⁿ(1,2),...,pⁿ(1,J),...,pⁿ(M,1),...,pⁿ(M,J)]^T表示系统状态在特定搜索深度n下的概率；系统状态转移概率Q＝q(m,j|m',j',Δt)；m,m'＝1,...,M；j,j'＝1,...,J由q(m,j|m',j',Δt)＝g(j|j',m',Δt)·h(m|m',j'→j,Δt)计算得到；其中，h(m|m',j'→j,Δt)表示部件状态配置在系统过程变量栅元j'→j转移过程中由m'变为m的概率；h(m|m',j'→j,Δt)为系统控制单元随机状态变化的全概率映射；假设部件失效之间相互独立，h(m|m',j'→j,Δt)能够简单地表示成每一单独控制单元失效概率的乘积；g(j|j',m',Δt)表示在给定部件状态配置m'下和时间步长Δt内，系统过程变量由栅元j'进入栅元j的条件概率；g(j|j',m',Δt)矩阵数据值能够通过系统控制数学模型或仿真模型/模拟器获得；通过对系统状态转移概率矩阵Q的计算，实现对系统状态概率映射关系模型的构建；

所述系统状态的始发点随机抽样方法为：通过一随机抽样函数每次从指定栅元V_j中随机抽取一个或多个数值作为系统状态的位置起始点用于系统轨迹模拟生成，随机抽样函数服从均匀分布；系统测试时，以栅元中的点代表栅元。

3.根据权利要求2所述的一种动态系统状态概率映射矩阵多向搜索方法，其特征在于：所述g(j|j',m',Δt)矩阵数据值能够通过系统控制数学模型采用连续积分法或仿真模型/模拟器采用等权点积分获得；此处的g(j|j',m',Δt)元素值通过仿真模型模拟得到。

4.根据权利要求1或2所述的一种动态系统状态概率映射矩阵多向搜索方法，其特征在于：系统仿真测试通过命令程序语句实现，测试结果包括系统落入栅元点的位置、数字化部件故障模式的符号参数以及最大和最小失效指数，并将所述测试结果存入数据库中。

5.根据权利要求2所述的一种动态系统状态概率映射矩阵多向搜索方法，其特征在于：所述系统状态条件转移概率通过以下等权点积分方法近似得到：假若系统测试过程中在栅元j'中随机抽取了u_j'个数作为系统状态起始位置点，这u_j'个数均匀分布在栅元j'中，经系统状态发展轨迹模拟，如果在时间步长Δt内有w_j个系统轨迹点落入栅元j，则系统状态条件转移概率矩阵近似于

6.根据权利要求2所述的一种动态系统状态概率映射矩阵多向搜索方法，其特征在于：所述多向矩阵动态搜索分析以系统指定初始状态s为起点，通过搜索系统状态转移概率矩阵Q中的非零元素，并将观测到的非零元素作为父节点继续下一步搜索，迭代循环直至指定搜索深度n；前、后搜索同步进行，每一搜索步都涉及符合元素的状态转移概率矩阵Q的更新计算，实行广度搜索优先，以遍历所有可能的序列路径；

其中，系统状态演变序列路径通过不同搜索步上系统状态栅元点前后相连构成的事件序列线路生成，系统状态演变序列路径的概率为每一搜索步分支概率的乘积；系统重要序列路径筛选方法如下：系统状态序列路径的重要性体现在其概率值排序上，通过截断准则ε以筛选出系统状态发展重要序列路径以及系统失效重要风险路径。

7.根据权利要求2所述的一种动态系统状态概率映射矩阵多向搜索方法，其特征在于：所述系统状态在不同时刻的可靠度计算方法如下：特定搜索深度n下，对系统状态进行分组归类，属于同一系统状态的序列分支归为一组，系统状态在n·Δt时刻的可靠度为组内所有序列分支概率之和，系统状态概率的预测通过Pⁿ⁺¹＝Q·Pⁿ实现。

Images