QNR Simulation

We illustrate here how to use the Variational Quantum Linear Solver as part of the Newton Raphson to simulate the water pressure and flow of a small water network using wntr-quantum.

Set up water network model

In this example, we test our quantum solvers into a slightly larger network as contained in Net2Loops.inp. Let's start by setting up the model:

import os
from pathlib import Path
import wntr
import wntr_quantum

# Define the directory path
epanet_tmp_dir = Path("/Users/murilo/scratch_dir/.epanet_quantum")

# check if the directory exists
if not epanet_tmp_dir.exists():
    epanet_tmp_dir.mkdir(parents=True, exist_ok=True)

# set the environment variables
os.environ["EPANET_TMP"] = str(epanet_tmp_dir)
os.environ["EPANET_QUANTUM"] = "/Users/murilo/Documents/NLeSC_Projects/Vitens/EPANET"

# set up network model
inp_file = "networks/Net2Loops_modified.inp"  # reproduces the values in https://doi.org/10.3390/w14060851
wn = wntr.network.WaterNetworkModel(inp_file)

# plot network
wntr.graphics.plot_network(wn, title=wn.name, node_labels=True)

# explicitly set Hazen-Williams head loss formulas and pressure dependent demand mode
wn.options.hydraulic.headloss = "H-W"  # 'D-W'
wn.options.hydraulic.demand_model = "PDD"  # 'PDD'

# required pressure: the pressure above which the consumer should receive the desired demand
wn.options.hydraulic.required_pressure = 30.0  # m
# minimum pressure: this is the pressure below which the consumer cannot receive any water
wn.options.hydraulic.minimum_pressure = 0.0  # m

# set time duration of the simulation (in seconds)
wn.options.time.duration = 0

# print options
dict(wn.options.hydraulic)
dict(wn.options.time)

import os from pathlib import Path import wntr import wntr_quantum # Define the directory path epanet_tmp_dir = Path("/Users/murilo/scratch_dir/.epanet_quantum") # check if the directory exists if not epanet_tmp_dir.exists(): epanet_tmp_dir.mkdir(parents=True, exist_ok=True) # set the environment variables os.environ["EPANET_TMP"] = str(epanet_tmp_dir) os.environ["EPANET_QUANTUM"] = "/Users/murilo/Documents/NLeSC_Projects/Vitens/EPANET" # set up network model inp_file = "networks/Net2Loops_modified.inp" # reproduces the values in https://doi.org/10.3390/w14060851 wn = wntr.network.WaterNetworkModel(inp_file) # plot network wntr.graphics.plot_network(wn, title=wn.name, node_labels=True) # explicitly set Hazen-Williams head loss formulas and pressure dependent demand mode wn.options.hydraulic.headloss = "H-W" # 'D-W' wn.options.hydraulic.demand_model = "PDD" # 'PDD' # required pressure: the pressure above which the consumer should receive the desired demand wn.options.hydraulic.required_pressure = 30.0 # m # minimum pressure: this is the pressure below which the consumer cannot receive any water wn.options.hydraulic.minimum_pressure = 0.0 # m # set time duration of the simulation (in seconds) wn.options.time.duration = 0 # print options dict(wn.options.hydraulic) dict(wn.options.time)

{'duration': 0.0,
 'hydraulic_timestep': 3600,
 'quality_timestep': 300,
 'rule_timestep': 360,
 'pattern_timestep': 7200,
 'pattern_start': 0.0,
 'report_timestep': 3600,
 'report_start': 0.0,
 'start_clocktime': 0.0,
 'statistic': 'NONE',
 'pattern_interpolation': False}

Solve model classically using EpanetSimulator

# define classical solver using EpanetSimulator
sim = wntr.sim.EpanetSimulator(wn)

# run the classical EPANET simulation
results_original_epanet = sim.run_sim()

# plot network
wntr.graphics.plot_network(
    wn,
    node_attribute=results_original_epanet.node["pressure"].iloc[0],
    link_attribute=results_original_epanet.link["flowrate"].iloc[0],
    node_colorbar_label="Pressure (m)",
    link_colorbar_label="Flow (liters/sec)",
    node_size=100,
    node_labels=True,
    link_labels=True,
    node_alpha=0.5,
)

results_original_epanet.node["pressure"], results_original_epanet.link["velocity"]

# define classical solver using EpanetSimulator sim = wntr.sim.EpanetSimulator(wn) # run the classical EPANET simulation results_original_epanet = sim.run_sim() # plot network wntr.graphics.plot_network( wn, node_attribute=results_original_epanet.node["pressure"].iloc[0], link_attribute=results_original_epanet.link["flowrate"].iloc[0], node_colorbar_label="Pressure (m)", link_colorbar_label="Flow (liters/sec)", node_size=100, node_labels=True, link_labels=True, node_alpha=0.5, ) results_original_epanet.node["pressure"], results_original_epanet.link["velocity"]

(name         2          3         4          5         6          7  \
 0     53.24765  30.466589  43.45023  33.808372  30.44595  30.553503   
 
 name             1  
 0     4.394531e-07  ,
 name         1         2         3         4         5        6         7  \
 0     1.894879  1.846526  1.462826  1.115641  1.136156  1.09944  1.298481   
 
 name         8  
 0     0.315217  )

Solve model using the classical Cholesky solver from QuantumEpanetSimulator

We now solve the same problem using the classical Epanet simulator. Note that, by default, QuantumEpanetSimulator uses a classical CholeskySolver by default to iteratively solve the linear problem.

# define classical solver using QuantumEpanetSimulator
sim = wntr_quantum.sim.QuantumEpanetSimulator(wn)

# run the EPANET quantum simulation using classical Cholesky
results_epanet = sim.run_sim()

# remember to set up EPANET Quantum environment variables!
epanet_path = os.environ["EPANET_QUANTUM"]
epanet_tmp = os.environ["EPANET_TMP"]

# check paths
print(f"Your EPANET quantum path: {epanet_path}")
print(f"Your EPANET temp dir: {epanet_tmp}\n")

# load EPANET A and b matrices from temp
epanet_A, epanet_b = wntr_quantum.sim.epanet.load_epanet_matrix()

# set the size of the Jacobian (A matrix)
epanet_A_dim = epanet_A.todense().shape[0]
print(f"Size of the Jacobian in EPANET simulator: {epanet_A_dim}")
print(f"Size of the b vector in EPANET simulator: {epanet_b.shape[0]}")

# save number of nodes and pipes
n_nodes = (len(results_epanet.node["pressure"].iloc[0]),)
n_pipes = len(results_epanet.link["flowrate"].iloc[0])

results_epanet.node["pressure"], results_epanet.link["velocity"]

# define classical solver using QuantumEpanetSimulator sim = wntr_quantum.sim.QuantumEpanetSimulator(wn) # run the EPANET quantum simulation using classical Cholesky results_epanet = sim.run_sim() # remember to set up EPANET Quantum environment variables! epanet_path = os.environ["EPANET_QUANTUM"] epanet_tmp = os.environ["EPANET_TMP"] # check paths print(f"Your EPANET quantum path: {epanet_path}") print(f"Your EPANET temp dir: {epanet_tmp}\n") # load EPANET A and b matrices from temp epanet_A, epanet_b = wntr_quantum.sim.epanet.load_epanet_matrix() # set the size of the Jacobian (A matrix) epanet_A_dim = epanet_A.todense().shape[0] print(f"Size of the Jacobian in EPANET simulator: {epanet_A_dim}") print(f"Size of the b vector in EPANET simulator: {epanet_b.shape[0]}") # save number of nodes and pipes n_nodes = (len(results_epanet.node["pressure"].iloc[0]),) n_pipes = len(results_epanet.link["flowrate"].iloc[0]) results_epanet.node["pressure"], results_epanet.link["velocity"]

Your EPANET quantum path: /Users/murilo/Documents/NLeSC_Projects/Vitens/EPANET
Your EPANET temp dir: /Users/murilo/scratch_dir/.epanet_quantum

Size of the Jacobian in EPANET simulator: 6
Size of the b vector in EPANET simulator: 6

(name          2          3          4          5         6          7  \
 0     53.247612  30.466719  43.450195  33.808483  30.44593  30.553486   
 
 name             1  
 0     4.394531e-07  ,
 name         1         2         3         4         5         6         7  \
 0     1.894883  1.846514  1.462829  1.115635  1.136152  1.099438  1.298483   
 
 name         8  
 0     0.315222  )

Check that the classical results are equivalent

import pandas.testing as pdt

pdt.assert_frame_equal(
    results_original_epanet.node["pressure"], results_epanet.node["pressure"]
)
pdt.assert_frame_equal(
    results_original_epanet.link["flowrate"], results_epanet.link["flowrate"]
)
# pdt.assert_frame_equal(results_original_epanet.link["velocity"], results_epanet.link["velocity"])

import pandas.testing as pdt pdt.assert_frame_equal( results_original_epanet.node["pressure"], results_epanet.node["pressure"] ) pdt.assert_frame_equal( results_original_epanet.link["flowrate"], results_epanet.link["flowrate"] ) # pdt.assert_frame_equal(results_original_epanet.link["velocity"], results_epanet.link["velocity"])

Solve water network with QuantumEpanetSimulator and VQLS

We now solve the model using VQLS. In this example, we are preconditioning the initial linear system using diagonal scaling and also using a mix of two classical optimizers.

import numpy as np

from qiskit.circuit.library import RealAmplitudes
from qiskit.primitives import Estimator
from qiskit_algorithms import optimizers as opt
from quantum_newton_raphson.vqls_solver import VQLS_SOLVER
from utils import compare_results

n_qubits = int(np.ceil(np.log2(epanet_A_dim)))

qc = RealAmplitudes(n_qubits, reps=3, entanglement="full")
estimator = Estimator()

linear_solver = VQLS_SOLVER(
    estimator=estimator,
    ansatz=qc,
    optimizer=[opt.COBYLA(maxiter=1000, disp=True), opt.CG(maxiter=500, disp=True)],
    matrix_decomposition="symmetric",
    verbose=True,
    preconditioner="diagonal_scaling",
    reorder=True,
)

sim = wntr_quantum.sim.QuantumEpanetSimulator(wn, linear_solver=linear_solver)
results_vqls = sim.run_sim(linear_solver=linear_solver)

results_classical, results_quantum = compare_results(results_epanet, results_vqls)

results_vqls.node["pressure"], results_vqls.link["velocity"]

import numpy as np from qiskit.circuit.library import RealAmplitudes from qiskit.primitives import Estimator from qiskit_algorithms import optimizers as opt from quantum_newton_raphson.vqls_solver import VQLS_SOLVER from utils import compare_results n_qubits = int(np.ceil(np.log2(epanet_A_dim))) qc = RealAmplitudes(n_qubits, reps=3, entanglement="full") estimator = Estimator() linear_solver = VQLS_SOLVER( estimator=estimator, ansatz=qc, optimizer=[opt.COBYLA(maxiter=1000, disp=True), opt.CG(maxiter=500, disp=True)], matrix_decomposition="symmetric", verbose=True, preconditioner="diagonal_scaling", reorder=True, ) sim = wntr_quantum.sim.QuantumEpanetSimulator(wn, linear_solver=linear_solver) results_vqls = sim.run_sim(linear_solver=linear_solver) results_classical, results_quantum = compare_results(results_epanet, results_vqls) results_vqls.node["pressure"], results_vqls.link["velocity"]

VQLS Iteration 1000 Cost 6.643e-03
   Return from subroutine COBYLA because the MAXFUN limit has been reached.

   NFVALS = 1000   F = 6.643069E-03    MAXCV = 0.000000E+00
   X =-6.003301E-01  -1.100831E+00  -1.251520E+00  -1.103777E+00   1.332641E-01
       2.956498E+00   8.709108E-01   1.550173E+00   7.461350E-01   2.731111E-01
      -2.069924E+00   1.162553E+00
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 74
         Function evaluations: 1599
         Gradient evaluations: 123
VQLS Iteration 1000 Cost 1.396e-02
   Return from subroutine COBYLA because the MAXFUN limit has been reached.

   NFVALS = 1000   F = 1.395057E-02    MAXCV = 0.000000E+00
   X = 4.137658E-01  -1.559580E+00   7.397351E-01   1.718300E+00   2.247472E+00
      -2.605908E+00  -5.912157E-01  -1.401045E+00   2.133813E+00   4.532001E+00
       1.376001E+00   2.235611E-01
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 60
         Function evaluations: 1456
         Gradient evaluations: 112
VQLS Iteration 1000 Cost 2.472e-03
   Return from subroutine COBYLA because the MAXFUN limit has been reached.

   NFVALS = 1000   F = 2.472077E-03    MAXCV = 0.000000E+00
   X = 3.857522E+00  -1.265118E+00   4.256407E+00  -1.239055E+00  -1.682544E+00
      -2.869380E+00  -4.414322E-01   4.111229E+00   3.386487E+00  -3.059907E+00
       2.906477E+00   2.793646E+00
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 51
         Function evaluations: 1196
         Gradient evaluations: 92
Quantum result 0.31111666560173035 within 0.008334552381341634% of classical result 0.31109073758125305
Quantum result 0.09354956448078156 within 0.016220675077726884% of classical result 0.09356474131345749
Quantum result 0.189796581864357 within 0.021948702387887445% of classical result 0.1897549331188202
Quantum result 0.009042594581842422 within 0.02503129426341868% of classical result 0.009044858627021313
Quantum result 0.14740082621574402 within 0.014701052975458082% of classical result 0.14737915992736816
Quantum result 0.05571400374174118 within 0.007790323802495752% of classical result 0.05570966377854347
Quantum result 0.06578942388296127 within 0.009217611617040038% of classical result 0.06579548865556717
Quantum result -0.00015969578817021102 within 0.019068349538275086% of classical result -0.000159726245328784
Quantum result 53.2465705871582 within 0.001955791657896612% of classical result 53.24761199951172
Quantum result 30.46950912475586 within 0.00915901406938471% of classical result 30.466718673706055
Quantum result 43.44719696044922 within 0.006900664149416603% of classical result 43.4501953125
Quantum result 33.812408447265625 within 0.011610469100186043% of classical result 33.8084831237793
Quantum result 30.44211769104004 within 0.012523151228295821% of classical result 30.44593048095703
Quantum result 30.54896354675293 within 0.014801334379934725% of classical result 30.553485870361328
Quantum result 4.39453117451194e-07 within 0.0% of classical result 4.39453117451194e-07

(name          2          3          4          5          6          7  \
 0     53.246571  30.469509  43.447197  33.812408  30.442118  30.548964   
 
 name             1  
 0     4.394531e-07  ,
 name         1         2        3         4         5         6         7  \
 0     1.895041  1.846214  1.46315  1.115356  1.136319  1.099524  1.298364   
 
 name         8  
 0     0.315162  )

import matplotlib.pyplot as plt

for i, result in enumerate(results_vqls.linear_solver_results):
    plt.semilogy(result.logger.values, label=f"iter # {i}")
    plt.legend()

import matplotlib.pyplot as plt for i, result in enumerate(results_vqls.linear_solver_results): plt.semilogy(result.logger.values, label=f"iter # {i}") plt.legend()

Plot Network model with absolute percent errors

To compare the results, we plot the network together with the absolute percent errors (with respect to the classical results) for the predicted pressures and flow rates.

from utils import get_ape_from_pd_series

wntr.graphics.plot_network(
    wn,
    node_attribute=get_ape_from_pd_series(
        results_vqls.node["pressure"].iloc[0], results_epanet.node["pressure"].iloc[0]
    ),
    link_attribute=get_ape_from_pd_series(
        results_vqls.link["flowrate"].iloc[0],
        results_epanet.link["flowrate"].iloc[0],
    ),
    node_colorbar_label="Pressure %",
    link_colorbar_label="Flows %",
    node_size=150,
    title=f"{inp_file}: Absolute Percent Error",
    node_labels=True,
    link_labels=True,
    node_alpha=0.5,
)

from utils import get_ape_from_pd_series wntr.graphics.plot_network( wn, node_attribute=get_ape_from_pd_series( results_vqls.node["pressure"].iloc[0], results_epanet.node["pressure"].iloc[0] ), link_attribute=get_ape_from_pd_series( results_vqls.link["flowrate"].iloc[0], results_epanet.link["flowrate"].iloc[0], ), node_colorbar_label="Pressure %", link_colorbar_label="Flows %", node_size=150, title=f"{inp_file}: Absolute Percent Error", node_labels=True, link_labels=True, node_alpha=0.5, )

<Axes: title={'center': 'networks/Net2Loops_modified.inp: Absolute Percent Error'}>

Plot pressures and flow rates

Let's check graphically the equivalence of the results.

import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

fig, ax = plt.subplots()
ax.scatter(
    results_classical[:n_pipes],
    results_quantum[:n_pipes],
    label="Flow rates",
    color="blue",
    marker="o",
)
ax.scatter(
    results_classical[n_pipes:],
    results_quantum[n_pipes:],
    label="Pressures",
    color="red",
    marker="s",
    facecolors="none",
)
ax.axline((0, 0), slope=1, linestyle="--", color="gray", label="")
ax.set_xlabel("Classical results")
ax.set_ylabel("Quantum results")
ax.legend()

# inset plot
ax_inset = inset_axes(ax, width="35%", height="35%", loc="lower right")
ax_inset.scatter(
    results_classical[:n_pipes], results_quantum[:n_pipes], color="blue", marker="o"
)
ax_inset.scatter(
    results_classical[n_pipes:],
    results_quantum[n_pipes:],
    color="red",
    marker="s",
    facecolors="none",
)
ax_inset.axline((0, 0), slope=1, linestyle="--", color="gray")
ax_inset.set_xlim(-0.25, 0.5)
ax_inset.set_ylim(-0.25, 0.5)
ax_inset.set_xticks([-0.25, 0, 0.5])
ax_inset.set_yticks([-0.25, 0, 0.5])
# plt.legend()
plt.show()

results_classical, results_quantum

import matplotlib.pyplot as plt from mpl_toolkits.axes_grid1.inset_locator import inset_axes fig, ax = plt.subplots() ax.scatter( results_classical[:n_pipes], results_quantum[:n_pipes], label="Flow rates", color="blue", marker="o", ) ax.scatter( results_classical[n_pipes:], results_quantum[n_pipes:], label="Pressures", color="red", marker="s", facecolors="none", ) ax.axline((0, 0), slope=1, linestyle="--", color="gray", label="") ax.set_xlabel("Classical results") ax.set_ylabel("Quantum results") ax.legend() # inset plot ax_inset = inset_axes(ax, width="35%", height="35%", loc="lower right") ax_inset.scatter( results_classical[:n_pipes], results_quantum[:n_pipes], color="blue", marker="o" ) ax_inset.scatter( results_classical[n_pipes:], results_quantum[n_pipes:], color="red", marker="s", facecolors="none", ) ax_inset.axline((0, 0), slope=1, linestyle="--", color="gray") ax_inset.set_xlim(-0.25, 0.5) ax_inset.set_ylim(-0.25, 0.5) ax_inset.set_xticks([-0.25, 0, 0.5]) ax_inset.set_yticks([-0.25, 0, 0.5]) # plt.legend() plt.show() results_classical, results_quantum

([0.31109074,
  0.09356474,
  0.18975493,
  0.009044859,
  0.14737916,
  0.055709664,
  0.06579549,
  -0.00015972625,
  53.247612,
  30.466719,
  43.450195,
  33.808483,
  30.44593,
  30.553486,
  4.3945312e-07],
 [0.31111667,
  0.093549564,
  0.18979658,
  0.009042595,
  0.14740083,
  0.055714004,
  0.065789424,
  -0.00015969579,
  53.24657,
  30.46951,
  43.447197,
  33.81241,
  30.442118,
  30.548964,
  4.3945312e-07])