【创新应用】新增 Lindblad 动力学的随机 Magnus 变分模拟算法

.gitignore

@@ -44,3 +44,4 @@ pyqpanda-algorithm/pyqpanda_alg/QSolver/*.lib
 /pyqpanda-algorithm/pyqpanda_alg/QSolver/libcurl.dll
 
 .ccls-cache
+test_outputs/

README.md

@@ -79,6 +79,9 @@ pyqpanda-algorithm 是由本源量子（Origin Quantum）开发的量子算法
 - **QSVD（量子变分奇异值分解）**  
   在变分框架下提取矩阵的奇异值与奇异向量，用于降维与推荐系统。
 
+- **LindbladMagnus（Lindblad 动力学的随机 Magnus 变分模拟）**  
+  基于随机 Magnus 展开（Scheme I–IV）与 McLachlan 变分原理的开量子系统变分模拟算法，内置 FMO 光合复合体、阻尼 TFIM、自由基对模型及 Liouvillian 精确解参考。
+
 ### 4. **通用工具与基础组件**
 
 提供量子振幅估计算法、比较器、稀疏编码等底层工具。

pyqpanda-algorithm/example/QAlgBase/testeg_Lindblad_FMO.py

@@ -0,0 +1,154 @@
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Demo: Lindblad dynamics of the Fenna-Matthews-Olson (FMO) pigment-protein
+complex, simulated through the stochastic Magnus expansion coupled to the
+McLachlan variational principle.
+
+The FMO complex is the canonical example of environment-assisted quantum
+transport: an electronic excitation is initialised on bacteriochlorophyll site
+1 and its population is transferred through the network of coupled sites
+towards a reaction centre (the sink) while being lost to the electronic ground
+state through radiative decay and dephased by the protein environment.
+
+The example follows the parameter set of
+    Huang et al., PRX Quantum 6, 040312 (2025),
+and contrasts the variational Lindblad-Magnus trajectory ensemble against the
+exact Liouvillian solution provided by :func:`pyqpanda_alg.LindbladMagnus.mesolve`.
+"""
+
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from pyqpanda_alg.LindbladMagnus import (HardwareEfficientAnsatz,
+                                         LindbladMagnusSolver, fmo_model,
+                                         mesolve)
+
+
+def main(out_dir: str = "test_outputs", traj_num: int = 30,
+         t_final: float = 200.0, n_steps: int = 41,
+         magnus_order: int = 1, layers: int = 2,
+         qsd_type: str = "nonlinear") -> None:
+    os.makedirs(out_dir, exist_ok=True)
+
+    # ------------------------------------------------------------------ #
+    #  Model: 5-site FMO sub-network (padded to 3 qubits)                #
+    # ------------------------------------------------------------------ #
+    H, c_ops, e_ops, psi0, labels = fmo_model()
+    print("=== Fenna-Matthews-Olson complex (5-site sub-network) ===")
+    print(f"Hamiltonian dimension: {H.shape[0]} (3 qubits, padded)")
+    print(f"Collapse operators:    {len(c_ops)} (3 dephasing + 3 decay + 1 sink)")
+    print(f"Observables:           {len(e_ops)} ({', '.join(labels)})")
+    print(f"Initial state:         |Site 1>")
+
+    # ------------------------------------------------------------------ #
+    #  Variational configuration                                         #
+    # ------------------------------------------------------------------ #
+    ansatz = HardwareEfficientAnsatz(n_qubits=3, layers=layers, init_state=psi0)
+    print(f"Ansatz:                HardwareEfficientAnsatz, "
+          f"{ansatz.n_parameters} parameters")
+
+    times = np.linspace(0.0, t_final, n_steps)
+    dt = times[1] - times[0]
+    print(f"Time grid:             {n_steps - 1} steps, "
+          f"dt = {dt:.3f}, T = {t_final:.1f}")
+
+    # ------------------------------------------------------------------ #
+    #  Exact reference (Liouvillian)                                     #
+    # ------------------------------------------------------------------ #
+    print("\nComputing exact Lindblad solution (Liouvillian) ...")
+    exact = mesolve(H, psi0, times, c_ops, e_ops)
+
+    # ------------------------------------------------------------------ #
+    #  Variational Lindblad-Magnus simulation                            #
+    # ------------------------------------------------------------------ #
+    print(f"\nRunning variational Lindblad-Magnus solver, "
+          f"order={magnus_order}, traj={traj_num}, {qsd_type} QSD ...")
+    solver = LindbladMagnusSolver(H, c_ops, ansatz,
+                                  qsd_type=qsd_type,
+                                  magnus_order=magnus_order,
+                                  integrator="rk4")
+    result = solver.solve(psi0, times, e_ops,
+                          traj_num=traj_num, seed=42, verbose=True)
+    mean = result.expect
+    std = result.std
+    max_err = float(np.abs(mean - exact).max())
+    print(f"  -> maximum error vs exact: {max_err:.4f}")
+    print(f"  -> final populations:")
+    for i, lab in enumerate(labels):
+        print(f"     {lab:>8}: exact = {exact[i, -1]:.4f}, "
+              f"sim = {mean[i, -1]:.4f} +- {std[i, -1]:.4f}")
+
+    # ------------------------------------------------------------------ #
+    #  Plotting                                                          #
+    # ------------------------------------------------------------------ #
+    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
+    palette = plt.cm.tab10.colors
+
+    ax = axes[0]
+    ax.set_xlabel("Time")
+    ax.set_ylabel("Population")
+    ax.set_xlim(times[0], times[-1])
+    for i, lab in enumerate(labels):
+        ax.plot(times, exact[i], "-", color=palette[i], label=f"Exact {lab}")
+        ax.plot(times, mean[i], "--", color=palette[i], label=f"Sim {lab}")
+    ax.set_title("FMO excitation transfer (Lindblad-Magnus variational)")
+    ax.legend(loc="best", fontsize=8)
+
+    ax = axes[1]
+    ax.set_xlabel("Time")
+    ax.set_ylabel("Absolute error")
+    ax.set_xlim(times[0], times[-1])
+    ax.set_yscale("log")
+    for i, lab in enumerate(labels):
+        err = np.abs(mean[i] - exact[i])
+        # Guard against log(0)
+        err_plot = np.where(err > 1e-12, err, 1e-12)
+        ax.plot(times, err_plot, color=palette[i], label=lab)
+    ax.set_title("Trajectory error vs exact Liouvillian")
+    ax.legend(loc="best", fontsize=8)
+
+    fig.tight_layout()
+    out_path = os.path.join(out_dir, "lindblad_magnus_fmo.png")
+    fig.savefig(out_path, dpi=200)
+    print(f"\nSaved figure to {out_path}")
+
+
+if __name__ == "__main__":
+    import argparse
+
+    parser = argparse.ArgumentParser(
+        description="FMO Lindblad-Magnus variational simulation demo")
+    parser.add_argument("--traj", type=int, default=30,
+                        help="Number of Lindblad trajectories (default 30).")
+    parser.add_argument("--t-final", type=float, default=200.0,
+                        help="Final simulation time (default 200).")
+    parser.add_argument("--n-steps", type=int, default=41,
+                        help="Number of time-grid points (default 41).")
+    parser.add_argument("--magnus-order", type=int, default=1,
+                        help="Order of the stochastic Magnus expansion "
+                             "(0-4, default 1).")
+    parser.add_argument("--layers", type=int, default=2,
+                        help="Number of hardware-efficient ansatz layers "
+                             "(default 2).")
+    parser.add_argument("--qsd-type", choices=["nonlinear", "linear"],
+                        default="nonlinear",
+                        help="QSD unravelling (default nonlinear).")
+    parser.add_argument("--out-dir", default="test_outputs",
+                        help="Output directory for figures.")
+    args = parser.parse_args()
+
+    main(out_dir=args.out_dir, traj_num=args.traj, t_final=args.t_final,
+         n_steps=args.n_steps, magnus_order=args.magnus_order,
+         layers=args.layers, qsd_type=args.qsd_type)

pyqpanda-algorithm/example/QAlgBase/testeg_Lindblad_TFIM.py

@@ -0,0 +1,108 @@
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Demo: Lindblad dynamics of the transverse-field Ising model with amplitude
+damping, simulated through the stochastic Magnus expansion coupled to the
+McLachlan variational principle.
+
+The example reproduces the qualitatively correct dynamics of the open TFIM
+system studied in
+    Huang et al., PRX Quantum 6, 040312 (2025),
+and compares the variational Lindblad-Magnus trajectories with the exact
+Liouvillian solution provided by :func:`pyqpanda_alg.LindbladMagnus.mesolve`.
+"""
+
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from pyqpanda_alg.LindbladMagnus import (HardwareEfficientAnsatz,
+                                         LindbladMagnusSolver, mesolve,
+                                         tfim_model)
+
+
+def main(out_dir: str = "test_outputs") -> None:
+    os.makedirs(out_dir, exist_ok=True)
+
+    # ------------------------------------------------------------------ #
+    #  Model: TFIM with two amplitude-damping channels                   #
+    # ------------------------------------------------------------------ #
+    H, c_ops, e_ops, psi0, labels = tfim_model()
+    print("=== Transverse-field Ising model with amplitude damping ===")
+    print(f"Hamiltonian dimension: {H.shape[0]} (2 qubits)")
+    print(f"Collapse operators:    {len(c_ops)}")
+    print(f"Observables:           {len(e_ops)} ({', '.join(labels)})")
+    print(f"Initial state:         |11>")
+
+    # ------------------------------------------------------------------ #
+    #  Variational configuration                                         #
+    # ------------------------------------------------------------------ #
+    # The hardware-efficient ansatz uses parameterised RX/RZ rotations on
+    # every qubit and parameterised RZZ entanglers in a ring.  All
+    # parameterised gates reduce to the identity at theta=0, so the ansatz
+    # reproduces the initial state |11> when theta = 0.
+    ansatz = HardwareEfficientAnsatz(n_qubits=2, layers=2, init_state=psi0)
+    print(f"Ansatz:                HardwareEfficientAnsatz, "
+          f"{ansatz.n_parameters} parameters")
+
+    # Time grid.  dt=0.05 is small enough to resolve the Ising dynamics
+    # (period ~ 2*pi/|g| ~ 6) and the moderate damping rate g=0.1.
+    times = np.linspace(0.0, 2.5, 51)
+    traj_num = 30
+
+    # ------------------------------------------------------------------ #
+    #  Exact reference (Liouvillian)                                     #
+    # ------------------------------------------------------------------ #
+    print("\nComputing exact Lindblad solution (Liouvillian) ...")
+    exact = mesolve(H, psi0, times, c_ops, e_ops)
+
+    # ------------------------------------------------------------------ #
+    #  Variational Lindblad-Magnus simulation                            #
+    # ------------------------------------------------------------------ #
+    fig_evo, ax_evo = plt.subplots(figsize=(8, 5))
+    ax_evo.set_xlabel("Time")
+    ax_evo.set_ylabel("Population")
+    ax_evo.set_xlim(times[0], times[-1])
+
+    palette = plt.cm.tab10.colors
+    for i, lab in enumerate(labels):
+        ax_evo.plot(times, exact[i], "-", color=palette[i],
+                    label=f"Exact {lab}")
+
+    for order in (1, 2):
+        print(f"\nRunning variational Lindblad-Magnus solver, "
+              f"order={order}, traj={traj_num} ...")
+        solver = LindbladMagnusSolver(H, c_ops, ansatz,
+                                      qsd_type="nonlinear",
+                                      magnus_order=order,
+                                      integrator="rk4")
+        result = solver.solve(psi0, times, e_ops,
+                              traj_num=traj_num, seed=42, verbose=True)
+        mean = result.expect
+        max_err = float(np.abs(mean - exact).max())
+        print(f"  -> maximum error vs exact: {max_err:.4f}")
+        for i, lab in enumerate(labels):
+            ax_evo.plot(times, mean[i], "--", color=palette[i],
+                        label=f"M{order} {lab}")
+
+    ax_evo.legend(loc="best", fontsize=8)
+    ax_evo.set_title("TFIM with amplitude damping: Lindblad-Magnus "
+                     "variational simulation")
+    out_path = os.path.join(out_dir, "lindblad_magnus_tfim.png")
+    fig_evo.tight_layout()
+    fig_evo.savefig(out_path, dpi=200)
+    print(f"\nSaved figure to {out_path}")
+
+
+if __name__ == "__main__":
+    main()

pyqpanda-algorithm/pyqpanda_alg/LindbladMagnus/README.md

@@ -0,0 +1,85 @@
+# LindbladMagnus — 开量子系统的变分量子模拟
+
+本模块实现了论文
+
+> J.-C. Huang, H.-E. Li, Y.-C. Wang, G.-Z. Zhang, J. Li, H.-S. Hu,  
+> *Towards Robust Variational Quantum Simulation of Lindblad Dynamics via
+> Stochastic Magnus Expansion*, **PRX Quantum** 6, 040312 (2025).
+> [arXiv:2503.22099](https://arxiv.org/abs/2503.22099)
+
+中的算法，将 **Lindblad 主方程** 通过量子态扩散（QSD）随机轨迹
+unravelling 后，用 **随机 Magnus 展开（Scheme I–IV）** 构造每一步的非厄米
+有效哈密顿量 `H_eff`，再用 **McLachlan 变分原理** 在参数化量子线路上完成
+一步变分时间演化。所有量子电路构造与态矢量模拟均基于 `pyqpanda3`，数值
+计算仅依赖 `numpy` / `scipy`，**不引入 qiskit / qutip 等其它量子框架**。
+
+## 快速上手
+
+```python
+import numpy as np
+from pyqpanda_alg.LindbladMagnus import (
+    LindbladMagnusSolver, HardwareEfficientAnsatz,
+    tfim_model, mesolve,
+)
+
+# 1. 准备开量子系统模型（这里用阻尼 TFIM 作为示例）
+H, c_ops, e_ops, psi0, labels = tfim_model()
+
+# 2. 构造变分 ansatz（theta=0 时为恒等，自动保留初态）
+ansatz = HardwareEfficientAnsatz(n_qubits=2, layers=2, init_state=psi0)
+
+# 3. 构造求解器并运行多条轨迹
+solver = LindbladMagnusSolver(H, c_ops, ansatz,
+                              qsd_type="nonlinear",
+                              magnus_order=1,     # 0=EM, 1..4=Magnus 阶数
+                              integrator="rk4")
+times = np.linspace(0, 2.5, 51)
+result = solver.solve(psi0, times, e_ops,
+                      traj_num=30, seed=42, verbose="tqdm")
+
+# 4. 与精确 Liouvillian 解对比
+exact = mesolve(H, psi0, times, c_ops, e_ops)
+print(f"最大误差: {np.abs(result.expect - exact).max():.4f}")
+```
+
+## 模块结构
+
+| 文件 | 职责 |
+| --- | --- |
+| `magnus.py` | 随机 Magnus 积分器（Scheme I–IV）与 Euler-Maruyama，生成 `H_eff` |
+| `ansatz.py` | `VariationalAnsatz` 与 `HardwareEfficientAnsatz`（参数化 RX/RY/RZ + RZZ/RXX/RYY） |
+| `variational.py` | 面向非厄米 `H_eff` 的 McLachlan 变分原理 + Euler / RK4 时间积分 |
+| `lindblad.py` | `LindbladMagnusSolver` 顶层求解器、`LindbladResult` 结果容器、函数式 `solve` |
+| `models.py` | FMO / TFIM / RPM 物理模型 + 自研 `mesolve` 精确解（Liouvillian） |
+
+## 关键 API
+
+- `LindbladMagnusSolver(...)` — 构造一次，多次调用 `.solve(...)` 复用配置
+- `LindbladResult` — dataclass，包含 `expect` / `std` / `norms` / `solver_info` / `seeds`
+- `solve(H, c_ops, ansatz, psi0, tlist, e_ops, ...)` — 一次性函数式接口
+- `mesolve(H, psi0, tlist, c_ops, e_ops)` — 基于 Liouvillian 的精确解（替代 `qutip.mesolve`）
+- `HardwareEfficientAnsatz(n_qubits, layers, init_state=...)` — 默认推荐 ansatz
+
+## 可选依赖
+
+- `tqdm` — 启用 `verbose="tqdm"` 进度条；未安装时自动回退到 `verbose=False`
+- `matplotlib` — 示例脚本绘图
+
+## 应用案例
+
+以下路径均相对仓库根目录：
+
+- `pyqpanda-algorithm/example/QAlgBase/testeg_Lindblad_TFIM.py` — TFIM 阻尼模型端到端 demo
+- `pyqpanda-algorithm/example/QAlgBase/testeg_Lindblad_FMO.py` — FMO 光合复合体能量传输 demo
+- `pyqpanda-algorithm/test/11-LindbladMagnus/demo01-LindbladMagnus-TFIM_FMO.ipynb` — notebook 教程
+- `test/LindbladMagnus/` — 单元测试与回归测试（`pytest`，已在 `test/pytest.ini` 注册）
+
+## 参数调优建议
+
+- **时间步长 `dt`**：建议从 `dt ≈ 0.05–0.1 × 2π/‖H‖` 起步；过大会导致 McLachlan
+  最小二乘方程病态，表现为 `Ground` / `Sink` 等通道人口过冲。
+- **Ansatz 层数 `layers`**：2–3 层对 TFIM 类问题足够；FMO 等多体系统建议 3+ 层
+  或换成 Hamiltonian Variational Ansatz。
+- **轨迹数 `traj_num`**：nonlinear QSD 下 20–50 条通常已收敛；linear QSD 下
+  因方差更大建议 ≥ 100 条。
+- **正则化 `eps`**：若日志中出现 `LinAlgError`，将 `eps` 调大到 `1e-8` 通常即可。