CN116706921B - 基于hhl算法的量子牛顿-拉夫逊法潮流计算方法和系统 - Google Patents

基于hhl算法的量子牛顿-拉夫逊法潮流计算方法和系统 Download PDF

Info

Publication number
CN116706921B
CN116706921B CN202310994633.8A CN202310994633A CN116706921B CN 116706921 B CN116706921 B CN 116706921B CN 202310994633 A CN202310994633 A CN 202310994633A CN 116706921 B CN116706921 B CN 116706921B
Authority
CN
China
Prior art keywords
node
voltage
quantum
power
ith
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202310994633.8A
Other languages
English (en)
Other versions
CN116706921A (zh
Inventor
张大波
徐震
马英浩
杨贺钧
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hefei University of Technology
Original Assignee
Hefei University of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hefei University of Technology filed Critical Hefei University of Technology
Priority to CN202310994633.8A priority Critical patent/CN116706921B/zh
Publication of CN116706921A publication Critical patent/CN116706921A/zh
Application granted granted Critical
Publication of CN116706921B publication Critical patent/CN116706921B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N10/00—Quantum computing, i.e. information processing based on quantum-mechanical phenomena
    • G06N10/60—Quantum algorithms, e.g. based on quantum optimisation, quantum Fourier or Hadamard transforms
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J3/00—Circuit arrangements for ac mains or ac distribution networks
    • H02J3/04—Circuit arrangements for ac mains or ac distribution networks for connecting networks of the same frequency but supplied from different sources
    • H02J3/06—Controlling transfer of power between connected networks; Controlling sharing of load between connected networks
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J2203/00—Indexing scheme relating to details of circuit arrangements for AC mains or AC distribution networks
    • H02J2203/20—Simulating, e g planning, reliability check, modelling or computer assisted design [CAD]
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Computing Systems (AREA)
  • Evolutionary Computation (AREA)
  • Condensed Matter Physics & Semiconductors (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Artificial Intelligence (AREA)
  • Power Engineering (AREA)
  • Supply And Distribution Of Alternating Current (AREA)

Abstract

本发明涉及电力系统潮流计算和量子计算技术领域，尤其是一种基于HHL算法的量子牛顿‑拉夫逊法潮流计算方法和系统。本发明提出的基于HHL算法的量子牛顿‑拉夫逊法潮流计算方法，结合量子计算和传统计算机两种计算资源进行电网潮流计算，本发明在有限量子硬件资源下构建的量子潮流计算方法，能够通过每次迭代生成新的修正方程降低HHL算法带来的误差，在不增加方程数和编程量的前提下，大大提高了现阶段量子交流潮流计算结果的精度。

Description

基于HHL算法的量子牛顿-拉夫逊法潮流计算方法和系统
技术领域
本发明涉及电力系统潮流计算和量子计算技术领域，尤其是一种基于HHL算法的量子牛顿-拉夫逊法潮流计算方法和系统。
背景技术
潮流计算是电力系统分析中最基本和关键的计算，能够依据电力系统网络中的参数和拓扑结构确定此系统的运行状态，被广泛应用于电力系统规划和调度运行中的各种研究。
现阶段复杂电力系统的潮流计算多采用快速解耦法、牛顿-拉夫逊法以及高斯-赛德尔迭代法。其中，在大规模电力系统潮流计算中，因直角坐标牛顿-拉夫逊法潮流计算的计算速度相比较慢，导致其使用受限。随着我国“双碳”目标的稳步推进，电网中的电气设备种类和数量陡然增长，电网的运行方式复杂多变，这使得直角坐标牛顿-拉夫逊法潮流计算的数据规模和计算时长陡然增长。针对当前模型驱动的潮流计算维度高并且实时性要求高的特点，传统计算方法和计算机的数据处理能力已接近极限，这成为制约构建新型电力系统的“卡脖子”问题。新型的计算方式和方法应被提出并应用于电网潮流计算。
量子计算作为一种新型的计算方式，利用量子叠加、量子纠缠等资源进行信息编码和处理，在海量信息存储和并行计算方面展现出其独有的优势，并已被证明在若干问题上具有相对于经典计算的极大优势。现阶段，量子计算已经成为内涵丰富的技术领域，涉及内容从最前沿的数学、物理等基础研究延伸到与诸多工程学科的交叉融合，再到高度工程化的应用技术开发，高速发展势头不减，并在数据快速搜索与排序、量子化学模拟、人工智能与机器学习等诸多领域表现出了可观的潜力，但在电力行业的应用仍处于起步阶段。
发明内容
为了克服上述现有技术中缺少通过量子计算实现电网潮流计算的技术的缺陷，本发明提出了一种基于HHL算法的量子牛顿-拉夫逊法潮流计算方法，利用量子计算机和传统计算机交替计算，提升了现有量子潮流计算的精度，即使面对现有有限的量子硬件，也可以保证较优的收敛次数和收敛精度。
本发明提出的一种基于HHL算法的量子牛顿-拉夫逊法潮流计算方法，包括以下步骤：
S1、结合电力系统的拓扑网络状态和运行参数，构建电力系统的节点导纳矩阵Y；统计电力系统中各节点有功功率P(i,0)和各节点无功功率Q(i,0)；设置各节点的电压，令电力系统拓扑网络中第i个节点的电压为V(i,1)，构建V(i,1)的复数表达形式作为V(i,1)表达式；i为序数，1≤i≤N，N为电力系统包含的平衡节点、PQ节点和PV节点的总数量；
S2、构建电力系统拓扑网络结构的非线性形式的量子牛顿-拉夫逊法潮流计算方程，简称非线性方程；
S3、将V(i,1)表达式和节点导纳矩阵Y代入非线性方程，通过方程求偏导，构建雅克比矩阵A，并构建不平衡量b；不平衡量b用于描述目标项的计算值和固定不变的已知值之间的偏差，目标项包括有功功率、无功功率和电压幅值中的一项或者多项；
S4、定义修正方程A×x=b，通过量子计算求解修正量x；
S5、根据修正量x更新非平衡节点的电压V(i,1)，并判断量子牛顿-拉夫逊法潮流计算是否达到收敛条件；是，则执行步骤S6；否，则判断t是否达到设定的迭代上限值；是，则执行步骤S6；否，则令t更新为t+1，再返回S3；
S6、输出所有节点的电压V(i,1)，完成量子牛顿-拉夫逊法潮流计算；其中，PQ节点和PV节点的电压V(i,1)为S5迭代更新后的电压V(i,1)，平衡节点电压V(i,1)为S1中设置的V(i,1)。
优选的，导纳矩阵Y中第i行第k列的元素Y(ik,1)为电力系统拓扑网络中第i个节点和第k个节点的互导纳，当i=k时，Y(ik,1)为第i个节点的自导纳；Y(ik,1)为复数；1≤k≤N；
Y(ik,1)=G(ik)+j×B(ik)
其中，G(ik)为Y(ik,1)的实部，B(ik)为Y(ik,1)的虚部，且G(ik)为电力系统拓扑网络中第i个节点和第k个节点之间电网支路的电导，B(ik)为电力系统拓扑网络中第i个节点和第k个节点之间电网支路的电纳；
非线性方程为：
V(i,1)×Σk∈[1,N][Y(ik,2)×V(k,2)]=P(i)+j×Q(i)
其中，V(k,1)表示电力系统拓扑网络中第k个节点的电压，V(k,2)为V(k,1)的共轭复数；Y(ik,2)为Y(ik,1)的共轭复数；P(i)为电力系统拓扑网络中第i个节点的有功功率计算值，Q(i)为电力系统拓扑网络中第i个节点的无功功率计算值，N为电力系统拓扑网络中的节点数量，Σk∈[1,N]表示求和时k的下限值为1且k的上限值为N；
S3中构建雅克比矩阵A的方式为：将V(i,1)的表达式和节点导纳矩阵Y代入非线性方程后，获取PQ节点和PV节点的有功功率P(i)的表示式以及PQ节点的无功功率Q(i)的表达式，令P(i)和Q(i)分别对V(i,1)表达式的自变量求偏导，根据偏导计算方程构建雅克比矩阵A。
优选的，S4中，首先将A转换为经典态A(1)，将b转换为经典态b(1)；A(1)为Hermitian矩阵，且A(1)为实对称矩阵，A(1)∈RM×M，M为设定值；b(1)为M维的单位向量，即b(1)=[b/|b|,0]T；然后重构修正方程并通过量子计算机求解重构的修正方程，以获取包含修正向量x的量子态x(1)，再通过量子逆编码获得x的解；
重构的修正方程为：
A(1)×x(1)=b(1)
b(1)=[b/|b|,0]T
x(1)=[0,x/|b|]T
其中，x(1)∈RM；AH为A的共轭转置矩阵。
优选的，采用直角坐标定义V(i,1)的表达式为：
V(i,1)=e(i)+j×f(i)
其中，e(i)为V(i,1)的实部，f(i)为V(i,1)的虚部；
P(i)的表达式为：
P(i)=Σk∈[1,N][e(i)×[G(ik)×e(k)-B(ik)×f(k)]+f(i)×[G(ik)×f(k)+B(ik)×e(k)]]
Q(i)的表达式为：
Q(i)=Σk∈[1,N][f(i)×[G(ik)×e(k)-B(ik)×f(k)]-e(i)×[G(ik)×f(k)+B(ik)×e(k)]]
其中，e(k)表示第k个节点的电压V(k,1)的实部，f(k)表示第k个节点的电压V(k,1)的虚部。
优选的，S3中构建的雅克比矩阵为：
其中，H、N、J、L、R和S均为雅克比矩阵的分块矩阵；
表示求偏导；H(ik)为分块矩阵H中第i行第k列的元素，N(ik)为分块矩阵N中第i行第k列的元素，J(ik)为分块矩阵J中第i行第k列的元素，L(ik)为分块矩阵L中第i行第k列的元素，R(ik)为分块矩阵R中第i行第k列的元素，S(ik)为分块矩阵S中第i行第k列的元素；
b=[ΔP，ΔQ，ΔV2]T
ΔP={ΔP(i)|i∈{PQ+PV}}T
ΔQ={ΔQ(i)|i∈{PQ}}T
ΔV2={ΔV(i)2|i∈{PV}}T
ΔP(i)=P(i,0)-Σk∈[1,N][e(i)×[G(ik)×e(k)-B(ik)×f(k)]+f(i)×[G(ik)×f(k)+B(ik)×e(k)]]
ΔQ(i)=Q(i,0)-Σk∈[1,N][f(i)×[G(ik)×e(k)-B(ik)×f(k)]-e(i)×[G(ik)×f(k)+B(ik)×e(k)]]
ΔV(i)2=V(i,0)2-V(i)2
P(i,0)为电力系统拓扑网络中第i个节点的有功功率，Q(i,0)为电力系统拓扑网络中第i个节点的无功功率，P(i,0)和Q(i,0)为电力系统已知的运行参数；ΔP(i)为第i个节点的有功功率计算值和已知值的误差，ΔP为非平衡节点的有功功率误差列向量，{PQ+PV}表示非平衡节点的集合；ΔQ(i)为第i个节点的无功功率的计算值和已知值的误差，ΔQ为PQ节点的无功功率误差列向量，{PQ}表示PQ节点的集合；V(i,0)为电力系统拓扑网络中第i个节点在S1中设置的电压幅值初始值，V(i)为电压V(i,1)的幅值，V(i)2=e(i)2+f(i)2；{PV}表示PV节点的集合；ΔV为PV节点的电压幅值误差向量，ΔV(i)为ΔV中第i个节点对应的误差；
x=[Δe，Δf]T
Δe={Δe(i)|i∈{PQ+PV}}T
Δf={Δf(i)|i∈{PQ+PV}}T
Δe表示电压实部的修正向量，Δe(i)为x中第i个节点的电压实部的修正值；Δf表示电压虚部的修正向量，Δf(i)为x中第i个节点的电压虚部的修正值。
优选的，S5中的收敛条件为：满足以下条件1a和2a中的至少一项；
1a：MAX||x|-|x’||<ε1；
x’为上一次迭代过程中求解的x；ε1为设定的精度；
x=[Δe，Δf]T
Δe={Δe(i)|i∈{PQ+PV}}T
Δf={Δf(i)|i∈{PQ+PV}}T
x’=[Δe’，Δf’]T
Δe’={Δe(i)’|i∈{PQ+PV}}T
Δf’={Δf(i)’|i∈{PQ+PV}}T
MAX||x|-|x’||=MAX{(||Δe(i)|-|Δe(i)’||，||Δf(i)|-|Δf(i)’||)|i∈{PQ+PV}}
Δe(i)’为x’中第i个节点的电压实部的修正值，Δf(i)’为x’中第i个节点的电压虚部的修正值；Δe’为x’中电压实部的修正向量，Δf’为x’中电压虚部的修正向量；
2a：MAX|x|<ε2；
MAX|x|=MAX{(|Δe(i)|，|Δf(i)|)|i∈{PQ+PV}}
ε2为设定的精度。
优选的，采用极坐标定义V(i,1)的表达式为：
V(i,1)=V(i)×ej×δ(i)=V(i)×[cosδ(i)+j×sinδ(i)]
其中，V(i)为V(i,1)的幅值，δ(i)为V(i,1)的相角；
P(i)的表达式为：
P(i)=V(i)×Σk∈[1,N][V(k)×[G(ik)×cosδ(ik)+B(ik)×sinδ(ik)]]
Q(i)的表达式为：
Q(i)=V(i)×Σk∈[1,N][V(k)×[G(ik)×sinδ(ik)-B(ik)×cosδ(ik)]]
其中，δ(ik)表示电力系统拓扑网络中第i个节点和第k个节点的相位差，即δ(ik)=δ(i)-δ(k)；δ(i)为V(i,1)的相角，δ(k)为V(k,1)的相角；V(k)表示第k个节点的电压V(k,1)的幅值。
优选的：
其中，H、N、J和L均为雅克比矩阵的分块矩阵；
其中，H(ik)为分块矩阵H中第i行第k列的元素，N(ik)为分块矩阵N中第i行第k列的元素，J(ik)为分块矩阵J中第i行第k列的元素，L(ik)为分块矩阵L中第i行第k列的元素；表示求偏导，V(k)表示第k个节点的电压V(k,1)的幅值，δ(k)表示第k个节点的电压V(k,1)的相角；
b=[ΔP，ΔQ]T
ΔP={ΔP(i)|i∈{PQ+PV}}T
ΔQ={ΔQ(i)|i∈{PQ}}T
ΔP(i)=P(i,0)-V(i)×Σk∈[1,N][V(k)×[G(ik)×cosδ(ik)+B(ik)×sinδ(ik)]]
ΔQ(i)=Q(i,0)-V(i)×Σk∈[1,N][V(k)×[G(ik)×sinδ(ik)-B(ik)×cosδ(ik)]]
其中，P(i,0)为电力系统拓扑网络中第i个节点的有功功率，Q(i,0)为电力系统拓扑网络中第i个节点的无功功率，P(i,0)和Q(i,0)为电力系统已知的运行参数；ΔP(i)为第i个节点的有功功率计算值和已知值的误差，ΔP为非平衡节点的有功功率误差列向量，{PQ+PV}表示非平衡节点的集合，即PQ节点和PV节点的集合；ΔQ(i)为第i个节点的无功功率的计算值和已知值的误差，ΔQ为PQ节点的无功功率误差列向量，{PQ}表示PQ节点的集合；
x=[Δδ，ΔV/V]T
Δδ={Δδ(i)|i∈{PQ+PV}}T
ΔV/V={ΔV(i)/V(i)|i∈{PQ}}T
Δδ表示电压相位的修正向量，Δδ(i)为x中第i个节点的电压相位的修正值；ΔV/V表示电压幅值的修正向量，ΔV(i)为第i个节点的电压幅值的修正值。
优选的，S5中的收敛条件为满足以下条件1b和2b中的至少一项；
1b：MAX||x|-|x’||<ε3；
x’为上一次迭代过程中求解的x；ε3为设定的精度；
x=[Δδ，ΔV/V]T
Δδ={Δδ(i)|i∈{PQ+PV}}T
ΔV/V={ΔV(i)/V(i)|i∈{PQ}}T
x’=[Δδ’，(ΔV/V)’]T
Δδ’={Δδ(i)’|i∈{PQ+PV}}T
(ΔV/V)’={ΔV(i)’/V(i)’|i∈{PQ}}T
Δδ(i)’为x’中第i个节点的电压相位的修正值，ΔV(i)’/V(i)’为x’中第i个节点的电压幅值的修正向量；Δδ’为x’中的电压相位的修正向量，(ΔV/V)’为x’中的电压幅值的修正向量；V(i)’为上一次迭代计算中V(i,1)的幅值，ΔV(i)’为上一次迭代计算中第i个节点的电压幅值修正值；
MAX||x|-|x’||=MAX{C1,C2}
C1=MAX{(||Δδ(i)|-|Δδ(i)’||)|i∈{PQ+PV}}
C2=MAX{(||ΔV(i)/V(i)|-|ΔV(i)’/V(i)’||)|i∈{PQ}}
C1和C2均为过渡项；
2b：MAX|b’|<ε4；
MAX|b’|=MAX{C3,C4}
C3={(|ΔP(i)’|)|i∈{PQ+PV}}
C4={(|ΔQ(i)’|)|i∈{PQ}}
C3和C4均为过渡项；ε4为设定的精度；
b’为根据求解后的x更新后的不平衡量b，b’的获取方式为：令非平衡节点的δ(i)更新为δ(i)+Δδ(i)，将PQ节点的V(i)更新为V(i)+ΔV(i)，根据更新后的各节点的电压幅值和电压相位计算不平衡量b记作b’；
令b’=[ΔP’，ΔQ’]T
ΔP’={ΔP(i)’|i∈{PQ+PV}}T
ΔQ’={ΔQ(i)’|i∈{PQ}}T
ΔP’为更新后的有功功率误差向量，ΔP(i)’为结合b’更新后的第i个节点的有功功率的计算值和已知值的误差；ΔQ’为更新后的无功功率误差向量，ΔQ(i)’为结合b’更新后的第i个节点的无功功率的计算值和已知值的误差。
本发明还提出了一种基于HHL算法的量子牛顿-拉夫逊法潮流计算系统，用于承载上述方法，系统包括存储器和处理器，存储器中存储有计算机程序，处理器连接存储器，处理器用于执行所述计算机程序，所述计算机程序被执行时用于实现所述的基于HHL算法的量子牛顿-拉夫逊法潮流计算方法。
本发明的优点在于：
（1）本发明提出的一种基于HHL算法的量子牛顿-拉夫逊法潮流计算方法，结合量子计算和传统计算机两种计算资源进行电网潮流计算，本发明在有限量子硬件资源下构建的量子潮流计算方法，能够通过每次迭代生成新的修正方程降低HHL算法带来的误差，在不增加方程数和编程量的前提下，大大提高了现阶段量子交流潮流计算结果的精度。
（2）本发明实现了与经典潮流计算方法相似的迭代过程，并相对于现有的量子潮流计算方法实现了更快的收敛速度和更高的精度。
（3）本发明对于不同规模的配电网和输电网，均能够在更少的迭代次数下实现量子潮流计算的稳定收敛。对于节点系统呈病态的电力系统网络，本发明能够通过线性化手段将潮流计算问题重构，降低矩阵的奇异值对于HHL算法求解误差的影响，进而实现量子潮流在此类系统的可靠计算和稳定收敛，且保证量子潮流计算结果的较高精度。
（4）本发明进一步给出了直角坐标和极坐标下的计算方式，便于在不同场景和需求下的应用，保证了方法的灵活运用和推广。
附图说明
图1为本发明的基于HHL算法的量子牛顿-拉夫逊法潮流计算方法的流程图；
图2为IEEE 4节点网络拓扑图；
图3为量子极坐标牛顿-拉夫逊潮流计算与传统极坐标牛顿-拉夫逊潮流计算的节点电压幅值求解结果；
图4为量子极坐标牛顿-拉夫逊潮流计算与传统极坐标牛顿-拉夫逊潮流计算的节点电压相角求解结果。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本实施方式中，分别结合直角坐标和极坐标，对本发明给出的图1所示的基于HHL算法的量子牛顿-拉夫逊法潮流计算方法做进一步阐述。
本实施方式中，基于HHL算法的量子直角坐标牛顿-拉夫逊法潮流计算方法具体包括以下步骤SA1-SA6。
SA1、令电力系统包含1个平衡节点、N1个PQ节点和N2个PV节点，电力系统的节点总数为N=1+N1+N2；根据电力系统拓扑网络和电力系统中的设备运行状况，构建节点导纳矩阵Y，导纳矩阵Y中第i行第k列的元素Y(ik,1)为电力系统拓扑网络中第i个节点和第k个节点的互导纳，当i=k时，Y(ik,1)为第i个节点的自导纳；Y(ik,1)为复数；1≤i≤N；1≤k≤N；
Y(ik,1)=G(ik)+j×B(ik) (1)
其中，G(ik)为Y(ik,1)的实部，B(ik)为Y(ik,1)的虚部，且G(ik)为电力系统拓扑网络中第i个节点和第k个节点之间电网支路的电导，B(ik)为电力系统拓扑网络中第i个节点和第k个节点之间电网支路的电纳；
统计电力系统中各节点有功功率P(i,0)和各节点无功功率Q(i,0)；
设置各节点的电压，令电力系统拓扑网络中第i个节点的电压为V(i,1)，构建V(i,1)的复数表达形式作为V(i,1)表达式；
V(i,1)=e(i)+j×f(i) (2a)
其中，e(i)为V(i,1)的实部，f(i)为V(i,1)的虚部。
SA2、构建电力系统拓扑网络结构的非线性形式的量子牛顿-拉夫逊法潮流计算方程，简称非线性方程。
非线性方程为：
V(i,1)×Σk∈[1,N][Y(ik,2)×V(k,2)]=P(i)+j×Q(i) (3)
其中，V(k,1)表示电力系统拓扑网络中第k个节点的电压，V(k,2)为V(k,1)的共轭复数；Y(ik,2)为Y(ik,1)的共轭复数；P(i)为电力系统拓扑网络中第i个节点的有功功率计算值，Q(i)为电力系统拓扑网络中第i个节点的无功功率计算值；N为电力系统拓扑网络中的节点数量，Σk∈[1,N]表示求和时k的下限值为1且k的上限值为N。
结合以下步骤SA3-SA6，对PQ节点的电压V(i,1)和PV节点的电压V(i,1)进行迭代，在迭代过程中PV节点的电压V(i,1)的幅值V(i)保持不变，平衡节点的电压V(i,1)保持不变。
SA3、在第t次迭代过程中，首先将第t-1次迭代后的所有节点的电压V(i,1)代入以上公式(2a)和(3)；第一次迭代时，所有节点的电压V(i,1)为设定的初始值；将公式(1)、(2a)代入公式(3)后，获取所有PQ节点和PV节点的有功功率P(i)的表示式如公式(3.1a)所示，获取所有PQ节点的无功功率Q(i)的表达式如公式(3.2a)所示，令P(i)和Q(i)分别对V(i,1)表达式的自变量e(i)和f(i)求偏导，根据偏导计算方程构建雅克比矩阵A；并构建不平衡量b。
P(i)=Σk∈[1,N][e(i)×[G(ik)×e(k)-B(ik)×f(k)]+f(i)×[G(ik)×f(k)+B(ik)×e(k)]] (3.1a)
Q(i)=Σk∈[1,N][f(i)×[G(ik)×e(k)-B(ik)×f(k)]-e(i)×[G(ik)×f(k)+B(ik)×e(k)]] (3.2a)
其中，e(k)表示第k个节点的电压V(k,1)的实部，f(k)表示第k个节点的电压V(k,1)的虚部。
针对PV节点，其电压幅值V(i)如公式(3.3a)所示；
e(i)2+f(i)2=V(i)2 (3.3a)
SA3中，首先将(3.1a)(3.2a)(3.3a)分别对e(i)和f(i)求偏导，将求得的偏导公式舍弃高次方项，然后构建雅克比矩阵A；
其中，H、N、J、L、R和S均为雅克比矩阵的分块矩阵；
表示求偏导；H(ik)为分块矩阵H中第i行第k列的元素，N(ik)为分块矩阵N中第i行第k列的元素，J(ik)为分块矩阵J中第i行第k列的元素，L(ik)为分块矩阵L中第i行第k列的元素，R(ik)为分块矩阵R中第i行第k列的元素，S(ik)为分块矩阵S中第i行第k列的元素。
不平衡量b指的是目标项的计算值和固定不变的已知值之间的偏差。在直角坐标牛顿-拉夫逊法潮流计算中，目标项包括：PQ节点的有功功率、PQ节点的无功功率、PV节点的有功功率和PV节点的电压幅值。
b=[ΔP，ΔQ，ΔV2]T
ΔP={ΔP(i)|i∈{PQ+PV}}T
ΔQ={ΔQ(i)|i∈{PQ}}T
ΔV2={ΔV(i)2|i∈{PV}}T
ΔP(i)=P(i,0)-Σk∈[1,N][e(i)×[G(ik)×e(k)-B(ik)×f(k)]+f(i)×[G(ik)×f(k)+B(ik)×e(k)]]
(4.1a)
ΔQ(i)=Q(i,0)-Σk∈[1,N][f(i)×[G(ik)×e(k)-B(ik)×f(k)]-e(i)×[G(ik)×f(k)+B(ik)×e(k)]]
(4.2a)
ΔV(i)2=V(i,0)2-V(i)2 (4.3a)
P(i,0)为电力系统拓扑网络中第i个节点的有功功率，Q(i,0)为电力系统拓扑网络中第i个节点的无功功率，P(i,0)和Q(i,0)为电力系统已知的运行参数；ΔP(i)为第i个节点的有功功率计算值和已知值的误差，ΔP为非平衡节点的有功功率误差列向量，{PQ+PV}表示非平衡节点的集合，即PQ节点和PV节点的集合；
ΔQ(i)为第i个节点的无功功率的计算值和已知值的误差，ΔQ为PQ节点的无功功率误差列向量，{PQ}表示PQ节点的集合；
V(i,0)为电力系统拓扑网络中第i个节点在S1中设置的电压幅值初始值；ΔV为PV节点的电压幅值误差向量，ΔV(i)为ΔV中第i个节点对应的误差；{PV}表示PV节点的集合。
SA4、定义修正方程Ax=b，通过量子计算结合以下步骤SA41-SA45求解修正向量x；
x=[Δe，Δf]T (5.1a)
Δe={Δe(i)|i∈{PQ+PV}}T
Δf={Δf(i)|i∈{PQ+PV}}T
Δe表示电压实部的修正向量，Δe(i)为x中第i个节点的电压实部的修正值；Δf表示电压虚部的修正向量，Δf(i)为x中第i个节点的电压虚部的修正值。
SA41、首先将A转换为经典态A(1)，将b转换为经典态b(1)；A(1)为Hermitian矩阵，且A(1)为实对称矩阵，A(1)∈RM×M，M=4(N-1)，N为电力系统拓扑网络中的节点数量；b(1)为M维的单位向量。
具体的，b(1)=[b/|b|,0]T，b(1)中前M/2维度的元素为b/|b|，后M/2维度的元素为0；
(5.2a)
其中，AH为A的共轭转置矩阵；
令x(1)为包含修正向量x的量子态。
具体的，x(1)=[0,x/|b|]T，x(1)中前M/2维度的元素为0，后M/2维度的元素为x/|b|，x(1)∈RM。
结合A(1)和b(1)重构修正方程为：A(1)×x(1)=b(1)
即： (5.3a)
SA42、在量子计算机中搭建包含辅助量子比特、量子傅里叶寄存器和输入寄存器的量子线路，将经典态A(1)采用哈密顿编码转换为量子态，将经典态b(1)采用振幅编码转换为量子态，将量子态的b(1)投影到量子态的A(1)的特征向量空间；
SA43、依次在量子线路中构建相位估计、控制旋转和逆相位估计的量子线路，用以处理处于量子态的A(1)和b(1)；
SA44、对辅助量子比特施加量子局部测量操作，输出辅助量子比特为|1>的量子态；
SA45、将SA44所得辅助量子比特为|1>的量子态逆编码以完成修正方程组Ax=b的量子求解，并将所得结果输出到传统计算机中进行后续计算，即后续步骤S5和S6在传统计算机中进行。
SA5、更新非平衡节点的电压V(i,1)，具体将e(i)更新为e(i)+Δe(i)，将f(i)更新为f(i)+Δf(i)，然后判断量子直角坐标牛顿-拉夫逊法潮流计算是否达到收敛条件；是，则执行步骤SA6；否，则判断t是否达到设定的迭代上限值；是，则执行步骤SA6；否，则令t更新为t+1，再返回SA3，以便将更新后的非平衡节点的电压实部e(i)和电压虚部f(i)返回SA3进行下一次更新。
步骤SA5中对列向量e和列向量f进行更新时，保持PV节点的电压幅值V(i,0)固定不变；更新范围为：
e={e(i)|i∈{PQ+PV}}T
f={f(i)|i∈{PQ+PV}}T
SA5中将更新后的列向量e和f代入公式(2a)，计算非平衡节点的电压V(i,1)。
具体的，本步骤SA5中的收敛条件为：满足以下条件1a和2a中的至少一项；
1a：MAX||x|-|x’||<ε1；
x’为上一次迭代过程中求解的x；ε1为设定的精度；
x=[Δe，Δf]T
Δe={Δe(i)|i∈{PQ+PV}}T
Δf={Δf(i)|i∈{PQ+PV}}T
x’=[Δe’，Δf’]T
Δe’={Δe(i)’|i∈{PQ+PV}}T
Δf’={Δf(i)’|i∈{PQ+PV}}T
MAX||x|-|x’||=MAX{(||Δe(i)|-|Δe(i)’||，||Δf(i)|-|Δf(i)’||)|i∈{PQ+PV}}=
MAX{MAX{(||Δe(i)|-|Δe(i)’||)i∈{PQ+PV}},MAX{(||Δf(i)|-|Δf(i)’||)|i∈{PQ+PV}}}
Δe(i)’为x’中第i个节点的电压实部的修正值，Δf(i)’为x’中第i个节点的电压虚部的修正值；Δe’为x’中电压实部的修正向量，Δf’为x’中电压虚部的修正向量；
2a：MAX|x|<ε2；
MAX|x|=MAX{(|Δe(i)|，|Δf(i)|)|i∈{PQ+PV}}
=MAX{MAX{(|Δe(i)|)|i∈{PQ+PV}}，MAX{(|Δf(i)|)|i∈{PQ+PV}}}
ε2为设定的精度。
SA6、统计所有节点的复数形式电压V(i,1)，其中非平衡节点的电压为SA5中迭代更新后的电压V(i,1)，平衡节点电压V(i,1)为S1中设置的V(i,1)，完成直角坐标量子牛顿-拉夫逊法潮流计算。
即SA6的输出为：
{V(i,1)|1≤i≤N}
V(i,1)=e(i)+j×f(i)
针对PQ节点和PV节点，e(i)和f(i)为S5迭代更新后输出的e(i)和f(i)；
针对平衡节点，e(i)和f(i)为设定的初始值。
本实施方式提出的另一种基于HHL算法的量子极坐标牛顿-拉夫逊法潮流计算方法，包括以下步骤SB1-SB6。
SB1、令电力系统包含1个平衡节点、N1个PQ节点和N2个PV节点，电力系统的节点总数为N=1+N1+N2；根据电力系统拓扑网络和电力系统中的设备运行状况，构建节点导纳矩阵Y，导纳矩阵Y中第i行第k列的元素Y(ik,1)为电力系统拓扑网络中第i个节点和第k个节点的互导纳，当i=k时，Y(ik,1)为第i个节点的自导纳；Y(ik,1)为复数；1≤i≤N；1≤k≤N；
Y(ik,1)=G(ik)+j×B(ik) (1)
其中，G(ik)为Y(ik,1)的实部，B(ik)为Y(ik,1)的虚部，且G(ik)为电力系统拓扑网络中第i个节点和第k个节点之间电网支路的电导，B(ik)为电力系统拓扑网络中第i个节点和第k个节点之间电网支路的电纳；
统计电力系统中各节点有功功率P(i,0)、各节点无功功率Q(i,0)；
设置各节点的电压，令电力系统拓扑网络中第i个节点的电压为V(i,1)，构建V(i,1)的复数表达形式作为V(i,1)表达式；
V(i,1)=V(i)×ej×δ(i)=V(i)×[cosδ(i)+j×sinδ(i)] (2b)
其中，V(i)为V(i,1)的幅值，δ(i)为V(i,1)的相角。
SB2、构建电力系统拓扑网络结构的非线性形式的量子牛顿-拉夫逊法潮流计算方程，简称非线性方程。
非线性方程为：
V(i,1)×Σk∈[1,N][Y(ik,2)×V(k,2)]=P(i)+j×Q(i) (3)
其中，V(k,1)表示电力系统拓扑网络中第k个节点的电压，V(k,2)为V(k,1)的共轭复数；Y(ik,2)为Y(ik,1)的共轭复数；P(i)为电力系统拓扑网络中第i个节点的有功功率计算值，Q(i)为电力系统拓扑网络中第i个节点的无功功率计算值；N为电力系统拓扑网络中的节点数量，Σk∈[1,N]表示求和时k的下限值为1且k的上限值为N。
结合以下步骤SB3-SB6，对PQ节点的电压V(i,1)和PV节点的电压V(i,1)进行迭代，在迭代过程中PV节点的电压V(i,1)的幅值V(i)保持不变，平衡节点的电压V(i,1)保持不变；即PV节点的电压V(i,1)的幅值V(i)始终为最初设置的初始值V(i,0)；平衡节点的电压V(i,1)的幅值V(i)始终为最初设置的初始值V(i,0)，平衡节点的电压V(i,1)的相角δ(i)始终为最初设置的初始值δ(i,0)。
SB3、在进行第t次迭代过程中，首先将第t-1次迭代后的所有节点的电压V(i,1)代入以上公式(2b)和(3)；第一次迭代时，所有节点的电压V(i,1)为设定的初始值；将公式(1)、(2b)代入公式(3)后，获取所有PQ节点和PV节点的有功功率P(i)的表示式如公式(3.1b)所示，获取所有PQ节点的无功功率Q(i)的表达式如公式(3.2b)所示，令P(i)和Q(i)分别对V(i,1)表达式的自变量V(i)和δ(i)求偏导，根据偏导计算方程构建雅克比矩阵A；并构建不平衡量b。
P(i)=V(i)×Σk∈[1,N][V(k)×[G(ik)×cosδ(ik)+B(ik)×sinδ(ik)]] (3.1b)
Q(i)=V(i)×Σk∈[1,N][V(k)×[G(ik)×sinδ(ik)-B(ik)×cosδ(ik)]] (3.2b)
δ(ik)表示电力系统拓扑网络中第i个节点和第k个节点的相位差，即δ(ik)=δ(i)-δ(k)；δ(i)为V(i,1)的相角，δ(k)为V(k,1)的相角。
SB3中，首先将(3.1b)(3.2b)分别对V(i)和δ(i)求偏导，将求得的偏导公式舍弃高次方项，然后构建雅克比矩阵A；
其中，H、N、J和L均为雅克比矩阵的分块矩阵；
H(ik)为分块矩阵H中第i行第k列的元素，N(ik)为分块矩阵N中第i行第k列的元素，J(ik)为分块矩阵J中第i行第k列的元素，L(ik)为分块矩阵L中第i行第k列的元素；表示求偏导，V(k)表示第k个节点的电压V(k,1)的幅值。
不平衡量b指的是目标项的计算值和固定不变的已知值之间的偏差。在极坐标牛顿-拉夫逊法潮流计算中，目标项包括：PQ节点的有功功率、PQ节点的无功功率和PV节点的有功功率。
b=[ΔP，ΔQ]T
ΔP={ΔP(i)|i∈{PQ+PV}}T
ΔQ={ΔQ(i)|i∈{PQ}}T
ΔP(i)=P(i,0)-V(i)×Σk∈[1,N][V(k)×[G(ik)×cosδ(ik)+B(ik)×sinδ(ik)]] (4.1b)
ΔQ(i)=Q(i,0)-V(i)×Σk∈[1,N][V(k)×[G(ik)×sinδ(ik)-B(ik)×cosδ(ik)]] (4.2b)
P(i,0)为电力系统拓扑网络中第i个节点的有功功率，Q(i,0)为电力系统拓扑网络中第i个节点的无功功率，P(i,0)和Q(i,0)为电力系统已知的运行参数；ΔP(i)为第i个节点的有功功率计算值和已知值的误差，ΔP为非平衡节点的有功功率误差列向量，{PQ+PV}表示非平衡节点的集合，即PQ节点和PV节点的集合；ΔQ(i)为第i个节点的无功功率的计算值和已知值的误差，ΔQ为PQ节点的无功功率误差列向量，{PQ}表示PQ节点的集合。
SB4、定义修正方程Ax=b，通过量子计算结合以下步骤SB41-SB45求解修正向量x；
x=[Δδ，ΔV/V]T (5.1b)
Δδ={Δδ(i)|i∈{PQ+PV}}T
ΔV/V={ΔV(i)/V(i)|i∈{PQ}}T
Δδ表示电压相位的修正向量，Δδ(i)为x中第i个节点的电压相位的修正值；ΔV/V表示电压幅值的修正向量，ΔV(i)为第i个节点的电压幅值的修正值，V(i)表示V(i,1)的幅值。
SB41、首先将A转换为经典态A(1)，将b转换为经典态b(1)；A(1)为Hermitian矩阵，且A(1)为实对称矩阵，A(1)∈RM×M，M=2(N-1+N1)，N为电力系统拓扑网络中的节点数量，N1为PQ节点数量；b(1)为M维的单位向量。
具体的，b(1)=[b/|b|,0]T，b(1)中前M/2维度的元素为b/|b|，后M/2维度的元素为0；
(5.2b)
其中，AH为A的共轭转置矩阵；
令x(1)为包含修正向量x的量子态。
具体的，x(1)=[0,x/|b|]T，x(1)中前M/2维度的元素为0，后M/2维度的元素为x/|b|，x(1)∈RM。
结合A(1)和b(1)重构修正方程为：A(1)×x(1)=b(1)
(5.3b)
SB42、在量子计算机中搭建包含辅助量子比特、量子傅里叶寄存器和输入寄存器的量子线路，将经典态A(1)采用哈密顿编码转换为量子态，将经典态b(1)采用振幅编码转换为量子态，将量子态的b(1)投影到量子态的A(1)的特征向量空间；
SB43、依次在量子线路中构建相位估计、控制旋转和逆相位估计的量子线路，用以处理处于量子态的A(1)和b(1)；
SB44、对辅助量子比特施加量子局部测量操作，输出辅助量子比特为|1>的量子态；
SB45、将SB44所得输出辅助量子比特为|1>的量子态逆编码以完成修正方程组Ax=b的量子求解，并将所得结果输出到传统计算机中进行后续计算。
SB5、更新非平衡节点的电压V(i,1)，具体将非平衡节点的δ(i)更新为δ(i)+Δδ(i)，将PQ节点的V(i)更新为V(i)+ΔV(i)，判断量子极坐标牛顿-拉夫逊法潮流计算是否达到收敛条件；是，则执行步骤SB6；否，则判断t是否达到设定的迭代上限值；是，则执行步骤SB6；否，则令t更新为t+1，再返回SB4，即将更新后的非平衡节点的电压V(i,1)返回SB4。
步骤SB5中更新非平衡节点的电压V(i,1)时，具体为更新列向量δ和列向量V，更新时保持PV节点的电压幅值V(i)固定不变；更新范围为：
δ={δ(i)|i∈{PQ+PV}}T
V={V(i)|i∈{PQ}}T
SB5中将更新后的列向量δ和V代入公式列向量(2b)，以更新非平衡节点的电压V(i,1)。
SB5中的收敛条件为满足以下条件1b和2b中的至少一项；
1b：MAX||x|-|x’||<ε3；
x’为上一次迭代过程中求解的x；ε3为设定的精度；
x=[Δδ，ΔV/V]T
Δδ={Δδ(i)|i∈{PQ+PV}}T
ΔV/V={ΔV(i)/V(i)|i∈{PQ}}T
x’=[Δδ’，(ΔV/V)’]T
Δδ’={Δδ(i)’|i∈{PQ+PV}}T
(ΔV/V)’={ΔV(i)’/V(i)’|i∈{PQ}}T
Δδ(i)’为x’中第i个节点的电压相位的修正值，ΔV(i)’/V(i)’为x’中第i个节点的电压幅值的修正向量；Δδ’为x’中的电压相位的修正向量，(ΔV/V)’为x’中的电压幅值的修正向量；V(i)’为上一次迭代计算中V(i,1)的幅值，ΔV(i)’上一次迭代计算中第i个节点的电压幅值修正值；
MAX||x|-|x’||=MAX{C1,C2}
C1=MAX{(||Δδ(i)|-|Δδ(i)’||)|i∈{PQ+PV}}
C2=MAX{(||ΔV(i)/V(i)|-|ΔV(i)’/V(i)’||)|i∈{PQ}}
C1和C2均为过渡项；
2b：MAX|b’|<ε4；
MAX|b’|=MAX{C3,C4}
C3={(|ΔP(i)’|)|i∈{PQ+PV}}
C4={(|ΔQ(i)’|)|i∈{PQ}}
C3和C4均为过渡项；ε4为设定的精度；
b’为根据求解后的x更新后的不平衡量b，b’的获取方式为：令非平衡节点的δ(i)更新为δ(i)+Δδ(i)，将PQ节点的V(i)更新为V(i)+ΔV(i)，将更新后的各节点的电压幅值和电压相位代入公式(4.1b)和(4.2b)，获得更新后的不平衡量b记作b’；
令b’=[ΔP’，ΔQ’]T
ΔP’={ΔP(i)’|i∈{PQ+PV}}T
ΔQ’={ΔQ(i)’|i∈{PQ}}T
ΔP’为更新后的有功功率误差向量，ΔP(i)’为结合b’更新后的第i个节点的有功功率的计算值和已知值的误差；ΔQ’为更新后的无功功率误差向量，ΔQ(i)’为结合b’更新后的第i个节点的无功功率的计算值和已知值的误差。
SB6、统计所有节点的复数形式电压V(i,1)，其中非平衡节点的电压为SB5中迭代更新后的电压V(i,1)，平衡节点电压V(i,1)为SB1中设置的V(i,1)，完成极坐标量子牛顿-拉夫逊法潮流计算。
以下结合具体实施例，对上述的基于HHL算法的量子极坐标牛顿-拉夫逊法潮流计算方法进行验证。
本实施例，基于包含4个节点的IEEE（标准节点系统）系统进行仿真实验。该IEEE系统中，节点1为PQ节点，节点2和节点3为PV节点，节点4为平衡节点，系统拓扑网络如图2所示。
本实施例中，分别采用经典极坐标牛顿-拉夫逊法潮流计算方法（简称经典计算方法）、现有的量子潮流计算方法（量子快速解耦潮流计算，简称现有量子计算方法）和本发明提出的基于HHL算法的量子极坐标牛顿-拉夫逊法潮流计算（简称本发明方法）针对本实施例提出的IEEE系统进行潮流计算。
本实施例中，现有量子计算方法在传统计算机上模拟量子计算完成，本发明方法的量子计算部分也在传统计算机上模拟量子计算完成。
本发明方法执行时，设置最大迭代次数itmax为100，收敛精度ε3=ε4=10-6；设置雅克比矩阵A和不平衡向量b的阶数均为4；本发明方法构建的量子线路图中，量子求解器包含一个辅助量子比特，量子求解器中的量子傅里叶寄存器和输入寄存器均由若干个量子比特组成；本发明方法在每次迭代计算中，依据各雅克比矩阵A的最大特征值，以及矩阵A(1)的阶数设置求解器中量子比特的数量。
采用本发明方法计算时，令节点1即PQ节点的电压幅值为VQ,bus-1，节点1的电压相角为δQ,bus-1，节点2的电压相角为δQ,bus-2，节点3的电压相角为δQ,bus-3；
采用经典计算方法时，节点1的电压幅值为VTr,bus-1，节点1的电压相角为δTr,bus-1，节点2的电压相角为δTr,bus-2，节点3的电压相角为δTr,bus-3；
本发明方法和经典计算方法的结算结果对比如图3、图4所示。结果表明本发明方法相对于经典计算方法，两者收敛趋势和迭代过程几乎一致。
现有量子计算方法的计算精度为10-5，考虑经典计算方法的计算精度高于现有量子计算方法，在此综合对比现有量子计算方法、本发明方法、经典计算方法在计算结果偏差、收敛迭代次数和计算精度三方面上的差异。
1）计算结果偏差。与经典计算方法的计算结果相比，现有量子计算方法和本发明方法的计算结果偏差如表1所示。依据表1的结果，本发明在各计算值上的误差均小于现有量子计算方法，证明了本发明方法实现了更优的计算精度。
表1 量子潮流计算方法的计算结果误差分析
表1中，以经典计算方法计算的参数（电压幅值和电压相角）作为实际值，以待评估方法（现有量子计算方法和本发明方法）计算的参数作为计算值；
误差计算公式为：
误差=（计算值/实际值）-1。
2）收敛迭代次数。对于现有量子计算方法，算法迭代计算8次时，图2所示系统的潮流计算收敛；对于本发明方法，结合图3、图4所知，算法迭代计算3次时，图2所示系统的潮流计算收敛，且与经典计算方法的迭代计算次数一致。依据上述的结果，可见本发明方法相对于现有的量子计算方法能够在更少的迭代次数下实现量子潮流计算的稳定收敛，且与其对应的经典潮流计算的收敛迭代次数一致，证明了本发明方法的有效性。由于本发明方法中部分计算可通过量子计算执行，故而结合量子计算后，本发明方法的计算速度相对于只能在传统计算机上执行的经典计算方法实现了很大的提高，极大的降低了潮流计算所需时间。
3）计算精度。现有量子潮流计算方法无法达到10-6的精度，本发明方法在图2所示系统实际执行过程中，在最后一次迭代计算中，不平衡量的最大绝对值MAX|b’|小于10-6，结果表明，相对于现有量子计算方法，本发明方法的收敛精度更高。
综上所述，通过与现有量子计算和经典计算方法的对比，证明了本发明方法的有效性，且本发明方法相对于现有量子计算，实现了更快的收敛速度和更高的精度。
当然，对于本领域技术人员而言，本发明不限于上述示范性实施例的细节，而还包括在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现的相同或类似结构。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。
此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。
本发明未详细描述的技术、形状、构造部分均为公知技术。

Claims (10)

1.一种基于HHL算法的量子牛顿-拉夫逊法潮流计算方法，其特征在于，包括以下步骤：
S1、结合电力系统的拓扑网络状态和运行参数，构建电力系统的节点导纳矩阵Y；统计电力系统中各节点有功功率P(i,0)和各节点无功功率Q(i,0)；设置各节点的电压，令电力系统拓扑网络中第i个节点的电压为V(i,1)，构建V(i,1)的复数表达形式作为V(i,1)表达式；i为序数，1≤i≤N，N为电力系统包含的平衡节点、PQ节点和PV节点的总数量；
S2、构建电力系统拓扑网络结构的非线性形式的量子牛顿-拉夫逊法潮流计算方程，简称非线性方程；
S3、将V(i,1)表达式和节点导纳矩阵Y代入非线性方程，通过方程求偏导，构建雅克比矩阵A，并构建不平衡量b；不平衡量b用于描述目标项的计算值和固定不变的已知值之间的偏差，目标项包括有功功率、无功功率和电压幅值中的一项或者多项；
S4、定义修正方程A×x=b，通过量子计算求解修正量x；
S5、根据修正量x更新非平衡节点的电压V(i,1)，并判断量子牛顿-拉夫逊法潮流计算是否达到收敛条件；是，则执行步骤S6；否，则判断t是否达到设定的迭代上限值；是，则执行步骤S6；否，则令t更新为t+1，再返回S3；
S6、输出所有节点的电压V(i,1)，完成量子牛顿-拉夫逊法潮流计算；其中，PQ节点和PV节点的电压V(i,1)为S5迭代更新后的电压V(i,1)，平衡节点电压V(i,1)为S1中设置的V(i,1)。
2.如权利要求1所述的基于HHL算法的量子牛顿-拉夫逊法潮流计算方法，其特征在于，导纳矩阵Y中第i行第k列的元素Y(ik,1)为电力系统拓扑网络中第i个节点和第k个节点的互导纳，当i=k时，Y(ik,1)为第i个节点的自导纳；Y(ik,1)为复数；1≤k≤N；
Y(ik,1)=G(ik)+j×B(ik)
其中，G(ik)为Y(ik,1)的实部，B(ik)为Y(ik,1)的虚部，且G(ik)为电力系统拓扑网络中第i个节点和第k个节点之间电网支路的电导，B(ik)为电力系统拓扑网络中第i个节点和第k个节点之间电网支路的电纳；
非线性方程为：
V(i,1)×Σk∈[1,N][Y(ik,2)×V(k,2)]=P(i)+j×Q(i)
其中，V(k,1)表示电力系统拓扑网络中第k个节点的电压，V(k,2)为V(k,1)的共轭复数；Y(ik,2)为Y(ik,1)的共轭复数；P(i)为电力系统拓扑网络中第i个节点的有功功率计算值，Q(i)为电力系统拓扑网络中第i个节点的无功功率计算值，N为电力系统拓扑网络中的节点数量，Σk∈[1,N]表示求和时k的下限值为1且k的上限值为N；
S3中构建雅克比矩阵A的方式为：将V(i,1)的表达式和节点导纳矩阵Y代入非线性方程后，获取PQ节点和PV节点的有功功率P(i)的表示式以及PQ节点的无功功率Q(i)的表达式，令P(i)和Q(i)分别对V(i,1)表达式的自变量求偏导，根据偏导计算方程构建雅克比矩阵A。
3.如权利要求2所述的基于HHL算法的量子牛顿-拉夫逊法潮流计算方法，其特征在于，S4中，首先将A转换为经典态A(1)，将b转换为经典态b(1)；A(1)为Hermitian矩阵，且A(1)为实对称矩阵，A(1)∈RM×M，M为设定值；b(1)为M维的单位向量；然后重构修正方程并通过量子计算机求解重构的修正方程，以获取包含修正向量x的量子态x(1)，再通过量子逆编码获得x的解；
重构的修正方程为：
A(1)×x(1)=b(1)
b(1)=[b/|b|,0]T
x(1)=[0,x/|b|]T
其中，x(1)∈RM；AH为A的共轭转置矩阵。
4.如权利要求3所述的基于HHL算法的量子牛顿-拉夫逊法潮流计算方法，其特征在于，采用直角坐标定义V(i,1)的表达式为：
V(i,1)=e(i)+j×f(i)
其中，e(i)为V(i,1)的实部，f(i)为V(i,1)的虚部；
P(i)的表达式为：
P(i)=Σk∈[1,N][e(i)×[G(ik)×e(k)-B(ik)×f(k)]+f(i)×[G(ik)×f(k)+B(ik)×e(k)]]
Q(i)的表达式为：
Q(i)=Σk∈[1,N][f(i)×[G(ik)×e(k)-B(ik)×f(k)]-e(i)×[G(ik)×f(k)+B(ik)×e(k)]]
其中，e(k)表示第k个节点的电压V(k,1)的实部，f(k)表示第k个节点的电压V(k,1)的虚部。
5.如权利要求4所述的基于HHL算法的量子牛顿-拉夫逊法潮流计算方法，其特征在于，S3中构建的雅克比矩阵为：
其中，H、N、J、L、R和S均为雅克比矩阵的分块矩阵；
；/>表示求偏导；H(ik)为分块矩阵H中第i行第k列的元素，N(ik)为分块矩阵N中第i行第k列的元素，J(ik)为分块矩阵J中第i行第k列的元素，L(ik)为分块矩阵L中第i行第k列的元素，R(ik)为分块矩阵R中第i行第k列的元素，S(ik)为分块矩阵S中第i行第k列的元素；
b=[ΔP，ΔQ，ΔV2]T
ΔP={ΔP(i)|i∈{PQ+PV}}T
ΔQ={ΔQ(i)|i∈{PQ}}T
ΔV2={ΔV(i)2|i∈{PV}}T
ΔP(i)=P(i,0)-Σk∈[1,N][e(i)×[G(ik)×e(k)-B(ik)×f(k)]+f(i)×[G(ik)×f(k)+B(ik)×e(k)]]
ΔQ(i)=Q(i,0)-Σk∈[1,N][f(i)×[G(ik)×e(k)-B(ik)×f(k)]-e(i)×[G(ik)×f(k)+B(ik)×e(k)]]
ΔV(i)2=V(i,0)2-V(i)2
P(i,0)为电力系统拓扑网络中第i个节点的有功功率，Q(i,0)为电力系统拓扑网络中第i个节点的无功功率，P(i,0)和Q(i,0)为电力系统已知的运行参数；ΔP(i)为第i个节点的有功功率计算值和已知值的误差，ΔP为非平衡节点的有功功率误差列向量，{PQ+PV}表示非平衡节点的集合；ΔQ(i)为第i个节点的无功功率的计算值和已知值的误差，ΔQ为PQ节点的无功功率误差列向量，{PQ}表示PQ节点的集合；V(i,0)为电力系统拓扑网络中第i个节点在S1中设置的电压幅值初始值，V(i)为电压V(i,1)的幅值，V(i)2=e(i)2+f(i)2；{PV}表示PV节点的集合；ΔV为PV节点的电压幅值误差向量，ΔV(i)为ΔV中第i个节点对应的误差；
x=[Δe，Δf]T
Δe={Δe(i)|i∈{PQ+PV}}T
Δf={Δf(i)|i∈{PQ+PV}}T
Δe表示电压实部的修正向量，Δe(i)为x中第i个节点的电压实部的修正值；Δf表示电压虚部的修正向量，Δf(i)为x中第i个节点的电压虚部的修正值。
6.如权利要求5所述的基于HHL算法的量子牛顿-拉夫逊法潮流计算方法，其特征在于，S5中的收敛条件为：满足以下条件1a和2a中的至少一项；
1a：MAX||x|-|x’||<ε1；
x’为上一次迭代过程中求解的x；ε1为设定的精度；
x=[Δe，Δf]T
Δe={Δe(i)|i∈{PQ+PV}}T
Δf={Δf(i)|i∈{PQ+PV}}T
x’=[Δe’，Δf’]T
Δe’={Δe(i)’|i∈{PQ+PV}}T
Δf’={Δf(i)’|i∈{PQ+PV}}T
MAX||x|-|x’||=MAX{(||Δe(i)|-|Δe(i)’||，||Δf(i)|-|Δf(i)’||)|i∈{PQ+PV}}
Δe(i)’为x’中第i个节点的电压实部的修正值，Δf(i)’为x’中第i个节点的电压虚部的修正值；Δe’为x’中电压实部的修正向量，Δf’为x’中电压虚部的修正向量；
2a：MAX|x|<ε2；
MAX|x|=MAX{(|Δe(i)|，|Δf(i)|)|i∈{PQ+PV}}
ε2为设定的精度。
7.如权利要求3所述的基于HHL算法的量子牛顿-拉夫逊法潮流计算方法，其特征在于，采用极坐标定义V(i,1)的表达式为：
V(i,1)=V(i)×ej×δ(i)=V(i)×[cosδ(i)+j×sinδ(i)]
其中，V(i)为V(i,1)的幅值，δ(i)为V(i,1)的相角；
P(i)的表达式为：
P(i)=V(i)×Σk∈[1,N][V(k)×[G(ik)×cosδ(ik)+B(ik)×sinδ(ik)]]
Q(i)的表达式为：
Q(i)=V(i)×Σk∈[1,N][V(k)×[G(ik)×sinδ(ik)-B(ik)×cosδ(ik)]]
其中，δ(ik)表示电力系统拓扑网络中第i个节点和第k个节点的相位差，即δ(ik)=δ(i)-δ(k)；δ(i)为V(i,1)的相角，δ(k)为V(k,1)的相角；V(k)表示第k个节点的电压V(k,1)的幅值。
8.如权利要求7所述的基于HHL算法的量子牛顿-拉夫逊法潮流计算方法，其特征在于：
其中，H、N、J和L均为雅克比矩阵的分块矩阵；
；其中，H(ik)为分块矩阵H中第i行第k列的元素，N(ik)为分块矩阵N中第i行第k列的元素，J(ik)为分块矩阵J中第i行第k列的元素，L(ik)为分块矩阵L中第i行第k列的元素；/>表示求偏导，V(k)表示第k个节点的电压V(k,1)的幅值，δ(k)表示第k个节点的电压V(k,1)的相角；
b=[ΔP，ΔQ]T
ΔP={ΔP(i)|i∈{PQ+PV}}T
ΔQ={ΔQ(i)|i∈{PQ}}T
ΔP(i)=P(i,0)-V(i)×Σk∈[1,N][V(k)×[G(ik)×cosδ(ik)+B(ik)×sinδ(ik)]]
ΔQ(i)=Q(i,0)-V(i)×Σk∈[1,N][V(k)×[G(ik)×sinδ(ik)-B(ik)×cosδ(ik)]]
其中，P(i,0)为电力系统拓扑网络中第i个节点的有功功率，Q(i,0)为电力系统拓扑网络中第i个节点的无功功率，P(i,0)和Q(i,0)为电力系统已知的运行参数；ΔP(i)为第i个节点的有功功率计算值和已知值的误差，ΔP为非平衡节点的有功功率误差列向量，{PQ+PV}表示非平衡节点的集合，即PQ节点和PV节点的集合；ΔQ(i)为第i个节点的无功功率的计算值和已知值的误差，ΔQ为PQ节点的无功功率误差列向量，{PQ}表示PQ节点的集合；
x=[Δδ，ΔV/V]T
Δδ={Δδ(i)|i∈{PQ+PV}}T
ΔV/V={ΔV(i)/V(i)|i∈{PQ}}T
Δδ表示电压相位的修正向量，Δδ(i)为x中第i个节点的电压相位的修正值；ΔV/V表示电压幅值的修正向量，ΔV(i)为第i个节点的电压幅值的修正值。
9.如权利要求8所述的基于HHL算法的量子牛顿-拉夫逊法潮流计算方法，其特征在于，S5中的收敛条件为满足以下条件1b和2b中的至少一项；
1b：MAX||x|-|x’||<ε3；
x’为上一次迭代过程中求解的x；ε3为设定的精度；
x=[Δδ，ΔV/V]T
Δδ={Δδ(i)|i∈{PQ+PV}}T
ΔV/V={ΔV(i)/V(i)|i∈{PQ}}T
x’=[Δδ’，(ΔV/V)’]T
Δδ’={Δδ(i)’|i∈{PQ+PV}}T
(ΔV/V)’={ΔV(i)’/V(i)’|i∈{PQ}}T
Δδ(i)’为x’中第i个节点的电压相位的修正值，ΔV(i)’/V(i)’为x’中第i个节点的电压幅值的修正向量；Δδ’为x’中的电压相位的修正向量，(ΔV/V)’为x’中的电压幅值的修正向量；V(i)’为上一次迭代计算中V(i,1)的幅值，ΔV(i)’为上一次迭代计算中第i个节点的电压幅值修正值；
MAX||x|-|x’||=MAX{C1,C2}
C1=MAX{(||Δδ(i)|-|Δδ(i)’||)|i∈{PQ+PV}}
C2=MAX{(||ΔV(i)/V(i)|-|ΔV(i)’/V(i)’||)|i∈{PQ}}
C1和C2均为过渡项；
2b：MAX|b’|<ε4；
MAX|b’|=MAX{C3,C4}
C3={(|ΔP(i)’|)|i∈{PQ+PV}}
C4={(|ΔQ(i)’|)|i∈{PQ}}
C3和C4均为过渡项；ε4为设定的精度；
b’为根据求解后的x更新后的不平衡量b，b’的获取方式为：令非平衡节点的δ(i)更新为δ(i)+Δδ(i)，将PQ节点的V(i)更新为V(i)+ΔV(i)，根据更新后的各节点的电压幅值和电压相位计算不平衡量b记作b’；
令b’=[ΔP’，ΔQ’]T
ΔP’={ΔP(i)’|i∈{PQ+PV}}T
ΔQ’={ΔQ(i)’|i∈{PQ}}T
ΔP’为更新后的有功功率误差向量，ΔP(i)’为结合b’更新后的第i个节点的有功功率的计算值和已知值的误差；ΔQ’为更新后的无功功率误差向量，ΔQ(i)’为结合b’更新后的第i个节点的无功功率的计算值和已知值的误差。
10.一种基于HHL算法的量子牛顿-拉夫逊法潮流计算系统，其特征在于，包括存储器和处理器，存储器中存储有计算机程序，处理器连接存储器，处理器用于执行所述计算机程序，所述计算机程序被执行时用于实现如权利要求1-9任一项所述的基于HHL算法的量子牛顿-拉夫逊法潮流计算方法。
CN202310994633.8A 2023-08-09 2023-08-09 基于hhl算法的量子牛顿-拉夫逊法潮流计算方法和系统 Active CN116706921B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310994633.8A CN116706921B (zh) 2023-08-09 2023-08-09 基于hhl算法的量子牛顿-拉夫逊法潮流计算方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310994633.8A CN116706921B (zh) 2023-08-09 2023-08-09 基于hhl算法的量子牛顿-拉夫逊法潮流计算方法和系统

Publications (2)

Publication Number Publication Date
CN116706921A CN116706921A (zh) 2023-09-05
CN116706921B true CN116706921B (zh) 2023-10-03

Family

ID=87829799

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310994633.8A Active CN116706921B (zh) 2023-08-09 2023-08-09 基于hhl算法的量子牛顿-拉夫逊法潮流计算方法和系统

Country Status (1)

Country Link
CN (1) CN116706921B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117035106B (zh) * 2023-10-10 2024-01-09 合肥工业大学 量子潮流计算修正方程的分阶计算方法、系统和存储介质
CN117977600B (zh) * 2024-04-02 2024-06-18 合肥工业大学 电力系统多潮流集成并行量子计算方法、系统和存储介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114268102A (zh) * 2021-12-23 2022-04-01 国网江苏省电力有限公司经济技术研究院 基于解析式概率潮流模型的电力系统运行状态量化方法
CN114759548A (zh) * 2022-04-02 2022-07-15 武汉大学 一种基于基追踪算法的配电网拓扑辨识与线损计算方法
CN115473233A (zh) * 2022-06-23 2022-12-13 深圳供电局有限公司 电压无功自律协同控制方法、装置和计算机设备
CN115954868A (zh) * 2022-12-19 2023-04-11 贵州电网有限责任公司 一种电网安全分析方法、装置及介质

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CA2827701A1 (en) * 2013-09-23 2015-03-23 Sureshchandra B. Patel Methods of patel decoupled loadlow computation for electrical power system

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114268102A (zh) * 2021-12-23 2022-04-01 国网江苏省电力有限公司经济技术研究院 基于解析式概率潮流模型的电力系统运行状态量化方法
CN114759548A (zh) * 2022-04-02 2022-07-15 武汉大学 一种基于基追踪算法的配电网拓扑辨识与线损计算方法
CN115473233A (zh) * 2022-06-23 2022-12-13 深圳供电局有限公司 电压无功自律协同控制方法、装置和计算机设备
CN115954868A (zh) * 2022-12-19 2023-04-11 贵州电网有限责任公司 一种电网安全分析方法、装置及介质

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
谢大为等."基于正负综合灵敏度的输电断面双层优化潮流控制策略".《电力工程技术》.2023,第42卷(第1期),全文. *

Also Published As

Publication number Publication date
CN116706921A (zh) 2023-09-05

Similar Documents

Publication Publication Date Title
CN116706921B (zh) 基于hhl算法的量子牛顿-拉夫逊法潮流计算方法和系统
CN113379057B (zh) 量子体系基态能量估计方法及系统
CN106356859B (zh) 一种基于Matlab的直角坐标牛顿法潮流计算方法
CN114580647B (zh) 量子系统的模拟方法、计算设备、装置及存储介质
US11175641B2 (en) Processing platform with holomorphic embedding functionality for power control and other applications
CN108804386B (zh) 一种电力系统负荷裕度的并行化计算方法
Eskandarpour et al. Experimental quantum computing to solve network DC power flow problem
CN107749627A (zh) 基于改进匹配追踪的智能配电网潮流雅可比矩阵估计方法
Bhutad et al. Three-phase load flow methods for radial distribution networks
Sun et al. Global state estimation for whole transmission and distribution networks
CN112600201A (zh) 基于多维全纯嵌入法的高维静态电压稳定边界计算方法
Leong et al. An LMI-based functional estimation scheme of large-scale time-delay systems with strong interconnections
Borzacchiello et al. Unified formulation of a family of iterative solvers for power systems analysis
CN108462181A (zh) 考虑稀疏性的智能配电网潮流雅可比矩阵鲁棒估计方法
CN117117842A (zh) 抗噪声的量子快速解耦潮流计算方法、系统和存储介质
CN111639463B (zh) 一种基于XGBoost算法的电力系统扰动后频率特征预测方法
CN114188945B (zh) 一种含光伏电源的配电网短路电流计算方法及装置
Ma et al. General optimization framework for accurate and efficient reconstruction of symmetric complex networks from dynamical data
CN111181166B (zh) 一种预测校正的不确定性仿射潮流方法
De Tuglie et al. Identification of dynamic voltage-current power system equivalents through artificial neural networks
CN109038588A (zh) 一种配电网潮流计算方法
Akhmetbaev et al. Investrgations of the Topological Method and Voltage Generation Algorithms for Nodes of Complex Electrical Networks
CN117977600B (zh) 电力系统多潮流集成并行量子计算方法、系统和存储介质
CN116739097B (zh) 量子测量设备性能估计方法及装置、电子设备和介质
CN109800540A (zh) 一种电力系统机电暂态仿真的初始化方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant