# Many-Speakers Single Channel Speech Separation
with Optimal Permutation Training

###### Abstract

Single channel speech separation has experienced great progress in the last
few years. However, training neural speech separation for a large number of
speakers (e.g., more than 10 speakers) is out of reach for the current
methods, which rely on the Permutation Invariant Training (PIT). In this work,
we present a permutation invariant training that employs the Hungarian
algorithm in order to train with an $O(C^{3})$ time complexity, where $C$ is
the number of speakers, in comparison to $O(C!)$ of PIT based methods.
Furthermore, we present a modified architecture that can handle the increased
number of speakers. Our approach separates up to $20$ speakers and improves
the previous results for large $C$ by a wide margin.

††*Equal contribution

Index Terms: speech separation, single channel, deep learning

##  1 Introduction

Single channel speech separation is a fundamental problem in the speech
process, which has seen tremendous advances in the last years. The main neural
architectures can be divided into two categories: (i) Spectral based [1, 2, 3]
and (ii) time domain based [4, 5, 6]. Currently, the latter category leads
with respect to the obtained accuracy.

Since the order of the speakers at the output of the neural network is
arbitrary, a permutation invariant training is performed. Most of the neural
architecture for speech separation use the permutation invariant training
(PIT) loss [7] or its extension to the utterance level (uPIT) [8]. Both
variants have a computational complexity of $O(C!)$, where $C$ is the number
of the speakers. As a result, it is not feasible to run PIT on more than ten
speakers.

In this work, following [9] and [10], we propose a novel method to train a
large number of speakers with a lower complexity of $O(C^{3})$, by using the
Hungarian algorithm. The Hungarian algorithm is able to find the optimal
permutation in terms of the minimal sum of pairwise losses, thus matching
between pairs of output- and target-signals. In order to enable the separation
network to deal with a large number of speakers, we further introduce an
architecture that combines two distinct approaches to separation networks,
LSTM and dilated convolutional layers.

In our experiments, our method separates up to $20$ speakers, which, as far as
we can ascertain, is twice the number tackled by any existing method.
Moreover, we show that our method improves the previous state of the art
separation results for separating $5$ and $10$ speakers.

##  2 Related Work

Signal channel speech separation was explored using classical approaches [11,
12] and, more recently, using deep learning methods. In [13] an LSTM neural
network with a phase sensitive loss function was introduced. An improvement in
SDR was demonstrated on the CHiME-2 [14] dataset. In [15] a neural separation
network with a clustering-based embedding was introduced, presenting results
for the separation of two speakers and introducing the WSJ-2mix dataset that
was extensively used by followup work. This work was further extended in [16]
by extracting an embedding of spectrogram segments and estimating a mask for
the separation part. Results were provided for two and three speakers and an
SDR improvement of $10.3$ dB and $7.1$ dB for WSJ-2mix and WSJ-3mix was
obtained. In [17], the neural separator network was introduced. Attractor
points in the embedding space were used to obtain the time-frequency bins for
each speaker. The improvement on the WSJ-mix dataset was by 5.49%.

Luo et al. [4] introduced TasNet, which is a time domain encoder-decoder
neural architecture for the single channel speech separation problem. They
show a results of 11.1 SDR improvement for WSJ-2mix dataset over the state of
the art. Wang et al. [3] proposed a neural architecture that separates the
speakers in both the time and the frequency domains simultaneously. They
presented an improvement of 13.2 SDR on the WSJ-2mix dataset. The work of [4]
further improved the architecture and introduced ConvTasNet [18], which
employed a dilated convolutional neural network, showing an SDR improvement
15.6 dB for the WSJ-2mix dataset.

Another improvement with LSTM network was introduced in [19], where the dual-
path recurrent neural network (DPRNN) architecture was employed to model
extremely long sequences. They showed SDR improvement of 18.08 dB on WSJ-2mix
dataset. In [5] a separation network with $MulCat$ blocks was introduced. The
proposed method also removed the masking sub-network, leading to an
improvement of 20.12 dB to WSJ-2mix dataset. Furthermore, the WSJ-mix dataset
was extended to include mixtures of $5$ speakers, where the SDR improvement
was 10.6 dB. [20] combined DPRNN and TasNet and for the WSJ-5mix dataset they
showed an SDR improvement of 10.41dB, and 11.14dB SDR improvement for online
remixing. Zeghidour et al. [6] introduced a neural separation network that
infers a representation to each speaker, by performing clustering, and used it
to separate the mixture. They show an SDR improvement of 22.2 dB for WSJ-2mix
dataset.

Since our work builds upon the $MulCat$ network architecture [5], we will
recap its major components. The network consists of an encoder, a separation
module and a decoder. The encoder and decoder are simple 1D convolutions. The
separation module starts with a chunking module which cuts the signal into
chunks in time. Then, a series of doubled $MulCat$ blocks is applied. During
training, a multi-scale loss is employed- after each doubled $MulCat$ block
the activations are reconstructed by the decoder into audio signals and fed
into the loss function. The method also uses the Scale-Invariant Signal-to-
Noise Ratio (SI-SNR) loss, which is a slight improvement to the traditional
SDR loss.

Another line of work for speech separation uses beamformers and introduces an
extension of the minimum variance distortionless response (MVDR) [21]. A
neural beamformer was introduced in [22] and further improved in [23] for the
speech separation problem. A follow up work introduced the linearly
constrained minimum variance (LCMV) beamformer [24].

In [25] a SinkPIT loss is introduced. They proposed a variant of the PIT loss,
which is based on Sinkhorn’s matrix balancing algorithm. They reduce the
complexity of the PIT loss from $O(C!)$ to $O(kC^{2})$, where $k$ is set to
$200$. It is important to note that the chosen permutation is only an
approximation of the optimal permutation. In another work, a probabilistic-PIT
loss which considers the output permutation as discrete latent random variable
was introduced [26]. [9] and [10] noted that it is possible to use the
Hungarian algorithm to find the best permutation for source separation.

###  2.1 Hungarian Algorithm

The linear sum assignment problem (also known as the assignment problem) is
the task of assigning $C$ agents to do $C$ tasks, such that each agent is
assigned to exactly one task, and the total cost for the agents performing the
tasks is minimal. In other words, given a $C$-by-$C$ matrix of costs, $M$, one
for each agent-task pair, find a permutation $\pi$ of the agents, such that
the sum of costs of paired agents and tasks is minimal:

$\pi=\underset{\pi\in\Pi_{C}}{\operatorname{argmin}}\sum_{i=1}^{n}M_{i,\pi(i)}$
(1)

A naive solution is to iterate over all $C!$ possible permutations.
Fortunately, an optimal and polynomial-time algorithm that solves the
assignment problem was proposed in $1955$ by Harold Kuhn [27], reviewed in
$1957$ by James Munkres [28] and is mostly known by the name the Hungarian
Algorithm. The initial time complexity of the algorithm was $O(C^{4})$ and it
was modified later on to a time complexity of $O(C^{3})$ [29].

Simply put, the algorithm starts by trying to find an obvious permutation. If
that fails, it goes on to make modifications to the input matrix, in order to
find a valid and optimal permutation. The number of the modification
iterations needed to find the solution is indicative of how well the
permutation fits the data compared to alternative permutations. In other
words, the more iterations needed until convergence, the less significant is
the optimal permutation.

##  3 Method

---
(i)
(ii)
Figure 1: (i) The proposed Hungarian Loss. $M$, a $C$-by-$C$ matrix of SI-SNR
losses between output and target pairs is computed. $M$ is fed into the
Hungarian algorithm, which efficiently finds the optimal permutation of target
signals, $\pi$. (ii) The proposed separation network architecture. The novel
components are the added Conv blocks and the Hungarian loss, which replaces
the PIT loss.

This work extends the work in [5], through a number of new contributions.
First, we introduce the Hungarian Loss which replaces the PIT loss and gives
an optimal solution to the permutation issue with a much lower time
complexity, $O(C^{3})$, which allows to train separation networks for many
speakers. Second, we introduce a new network architecture that uses stacked
dilated convolutions before each pair of $MulCat$ blocks of [5]. The overall
method, is depicted in Figure 1.

###  3.1 Hungarian Loss

A single-channel speech separation network takes an audio signal that contains
a mixture of $C$ speakers speaking and outputs $C$ audio signals, each
optimized to contain a separate speaker.

During training, the network outputs the separated audio signals in an
arbitrary order. Thus, in order to compute a meaningful loss, an alignment,
i.e. a permutation, needs to be recovered between the outputs of the network
and the separated target signals. One way to find the right permutation is to
iterate over all possible $C!$ permutations and choose the one which gives the
lowest mean loss value on the pairs (i.e. PIT). The computational cost of PIT
is unnoticeable when $C$ is small, in comparison to the other parts of the
network. However, it makes training on a large number of speakers impossible.
For instance, for 20 speakers PIT needs to check $20!\approx 2.4\times
10^{18}$ different permutations).

To address this issue, we formulate the task of finding the permutation which
minimizes the loss function as a linear sum assignment problem. Given the $C$
output signals and $C$ target signals, we calculate the pairwise loss value,
$\hat{\ell}(s_{i},\hat{s}_{j})$, on every pair of output ($\hat{s}_{j}$) and
target ($s_{i}$) signals, which gives an $C$-by-$C$ matrix of losses, $M$.
Next, we assign each output with a unique target and vice-versa. Such
assignment is equivalent to choosing $C$ elements of the matrix, such that
each chosen element is in a different row and column from all others. An
optimal assignment minimizes the sum of values of the chosen elements. The PIT
loss can then be viewed as a brute-force solution to this problem, iterating
over all possible solutions:

$\ell(s,\hat{s})=\min_{\pi\in\Pi_{C}}~{}\frac{1}{C}\sum_{i=1}^{C}\hat{\ell}(s_{i},\hat{s}_{\pi(i)})$
(2)

By running the Hungarian Algorithm on $M$, we efficiently find the optimal
permutation in polynomial time instead of the brute-force, factorial-time PIT.
Note that the assignment algorithm does not need to be differential, since we
find a permutation of the targets by which we calculate the loss, meaning that
the process is separate from the backwards calculation of gradients.

###  3.2 Model

We shift our focus to solving the task of separating mixtures of many ($C\geq
10$) speakers. As $C$ increases, the task of separating the mixtures becomes
more challenging. Thus, we propose a new and suitable network architecture.

For this end, we modify the $MulCat$-based architecture [5] by adding stacked
dilated convolutions before each pair of $MulCat$ blocks. In addition, we
increased some of the network’s hyper-parameters to achieve a larger capacity
needed for the harder tasks. The dilated convolutions scheme is borrowed from
[18]: We use $8$ 1-D Conv blocks stacked on top of each other, with dilation
factors $2^{i-1}$ for $i\in\\{1,2,...,8\\}$. This corresponds to a single
column in the separation module of [18].

The rest of the model is in accordance with [5], i.e. using the same encoder,
chunking, $MulCat$ blocks, decoder and SI-SNR loss. We also adopt the multi-
scale loss scheme, which applies the loss function after each double MulCat
block, instead of just the last. We tuned some hyper-parameters of the
architecture to accompany the harder tasks when $C$ is large: $N$, number of
features, was increased from $128$ to $256$. $L$, the encoder’s kernel size
was increased from $8$ to $16$. $H$, the number of hidden units in the LSTMs
was increased from $128$ to $256$ and finally $R$, the number of double
$MulCat$ blocks was increased from $6$ to $7$. A similar hyper-parameters
adjustment was done in [25].

Figure 2: The average number of iterations per example performed by the
Hungarian algorithm. The average is decreasing as training progresses, which
indicates that the separation is improving. This also means that the runtime
of the Hungarian algorithm is getting even shorter.
---
Figure 3: Separation results for a mixture of $10$ speakers. First row: mel spectrogram of ground truth signals. Second row: mel spectrogram of SinkPIT signals. Third row: mel spectrogram of the outputs of our method. The last row shows the SI-SDR improvement for the SinkPIT method and our method. The x-axis is sorted by the SI-SDRi in descending order. |
---|---
(a) | (b)
Figure 4: The pairwise SI-SDR matrix $M$, sorted in descending order. For the
same input sample as in Figure 3. (a) Our method. (b) SinkPIT [25].

##  4 Experiments

Comparison to state of the art  We show results on datasets derived from WSJ
corpus [30] and LibriSpeech [31]. For WSJ, we use the 5-speaker mix,
introduced in [5], which uses the same procedure as in [15], i.e. $30$ hours
of speech from the training set si$\\_$tr$\\_$s were used to create the
training and validation sets. The five speakers were randomly chosen and
combined with random SNR values between $0-5$ dB. The test set is created from
si$\\_$et$\\_$s and si$\\_$dt$\\_$s with 16 speakers, that differ from the
speakers of the training set. For LibriSpeech, we use the LibriMix [32]
datasets with mixes of 5, 10, 15 and 20 speakers. LibriMix offers mixtures of
2 and 3 speakers from LibriSpeech, and we used the given scripts to create
mixtures of 5, 10, 15 and 20 speakers. Some modifications to the script needed
to be done. These can be found online*. We used LibriMix’s given parameters to
get a sample rate of $8$Khz, clean mixtures (no noise added) and each sample
is cut in length, according to the minimal length sample in the mixture. We
also use the augmentation process as in [25].

††*https://github.com/ShakedDovrat/LibriMix

A separate model is trained for each dataset, with the corresponding number of
output channels. Training was done using the Adam optimizer [33], with batch
size $32$ and a learning rate of $1e-3$ which was multiplied by $0.95$ every
two epochs. During training, each sample is cut into 4-second segments. Table
1 compares the results of our model to other models, using the SI-SDRi
metrics. As can be seen, our model outperforms the other methods by a large
margin. For Libri-5mix and Libri-10mix, we improve the previous results by
$1.89$dB and $1.33$dB respectively. Interestingly, the previous state of the
art results for $Libri5Mix$ is obtained with $MulCut$ [5], whereas for
$Libri10Mix$ it is obtained by [25], which uses SinkPIT. This is due to the
fact that the PIT loss in $MulCut$ for 10 speakers is prohibitive. For
WSJ-5mix, the SDR of our method improves by more than $2$dB over the previous
method. We are the first to present results for the $Libri15Mix$ and
$Libri20Mix$ datasets and running previous work on it is not practical. Sample
results are shared online https://shakeddovrat.github.io/hungarian/.

Hungarian loss vs. alternative losses  In order to show the benefits of using
the Hungarian algorithm as opposed to PIT, we show the training duration of an
epoch of each dataset, depicted in Table 2. As shown, on a $5$ speaker mix
both methods take about the same time. However, for $C=10$, the Hungarian
method is about $9$ times faster on our model. For $C\geq 15$, the Hungarian
method is still fast, while PIT is unusable as it failed to complete a single
epoch after days of running. To our knowledge, the SinkPIT can run on $20$
speakers but only find an approximation to the optimal permutation, whereas
our method is both optimal and fast (the source code was not published).

In Figure 2, we present the average number of iterations per example that were
performed by the Hungarian algorithm for the various datasets. We can observe
two phenomena: (i) As the training proceeds and the neural network improves
the separation performance, the average number of iterations is decreased.
(ii) The average number of iterations is lower when the number of speakers in
the mixture is lower. Both of these observations came from the fact that the
Hungarian algorithm needs no iterations to converge when the outputs or the
network are mixed or noisy.

In Figure 3, we plot the mel spectrogram for typical samples as in [25]. This
is shown for a mixture of $10$ speakers, from the $Libri10Mix$ dataset. As can
be seen in speakers $4$ and $8$ our method provide a cleaner mel spectrogram,
Furthermore, the SI-SDRi difference is 4.9 dB and 2.4 dB for speakers $4$ and
$8$ respectively (there is a significant improvement in all speakers).

In Figure 4, we plot the pairwise SI-SDR negative matrix $M$ sorted in
descending order, in comparison to the results of the SinkPIT system [25].
Evidently, the entropy in our method is lower, especially in the last rows,
i.e., it has much less confusion compared to the baseline method.

Ablation study  We run ablation analysis in order to understand the
contribution of each component of our method. The results are summarized in
Table 3, where -C means without the 1D convolutions and -H means using the PIT
loss instead of the Hungarian loss. As can be seen, adding the dilated
convolution to the architecture improves the performance. Moreover, without
the Hungarian algorithm, it is practically impossible to train on $15$ or more
speakers. For $10$ speakers, using PIT drastically diminishes performance due
to longer training time. On the other hand, for $5$ speakers, the PIT and
Hungarian have similar performance since they both found the optimal
permutation in a similar runtime.

##  5 Conclusions

In this work, we provide a method for single channel sound separation for a
large number of sources. Our method is the first work to show that one can
separate a mixture of $20$ speakers from a single channel recording. Our
solution is based on the Hungarian algorithm, which efficiently finds the
optimal permutation from the $C!$ possible permutations and on a new network
architecture that adds stacked 1-D convolutions and added capacity to the
state of the art architecture.

Table 1: SDR improvement performance of various models versus number of speakers and datasets. ’X’ results indicated simulation that failed to complete a single epoch, due to high complexity PIT loss. xMix is LibriMix with x speakers. Model | 5Mix | 10Mix | 15Mix | 20Mix | WSJ-5
---|---|---|---|---|---
ConvTasNet | - | - | - | - | 6.8
DPRNN [19] | - | - | - | - | 8.6
MulCat [5] | 10.83 | 4.74 | X | X | 10.6
TasTas [20] | - | - | - | - | 11.14
SinkPIT [25] | 9.39 | 6.45 | - | - | -
Ours | 12.72 | 7.78 | 5.66 | 4.26 | 13.22
Table 2: Training run times of PIT and Hungarian algorithm in minutes per epoch. ’X’ indicates simulation that failed to complete a single epoch, due to long run times of the PIT loss. Dataset | #Spkrs | #Perms | PIT | Hungarian
---|---|---|---|---
WSJ-5mix | 5 | 120 | 65 | 62
Libri-5Mix | 5 | 120 | 139 | 140
Libri-10Mix | 10 | $\approx 3.6e6$ | 462 | 52
Libri-15Mix | 15 | $\approx 1.3e12$ | X | 36
Libri-20Mix | 20 | $\approx 2.4e18$ | X | 29
Table 3: Ablation analysis - SDR Performance for LibriMix and WSJ-mix datasets. xMix is LibriMix with x speakers. ”-C” means without the 1D convolutions and ”-H” means using PIT instead of the Hungarian loss. ’X’ indicates simulation that failed to complete a single epoch, due to long run times of the PIT loss. Model | 5Mix | 10Mix | 15Mix | 20Mix | WSJ-5
---|---|---|---|---|---
Ours-C-H | 11.10 | 4.47 | X | X | 12.18
Ours-H | 12.53 | 4.84 | X | X | 13.07
Ours-C | 11.19 | 5.89 | 5.14 | 4.10 | 12.10
Ours | 12.72 | 7.78 | 5.66 | 4.26 | 13.22

##  6 Acknowledgments

This project has received funding from the European Research Council (ERC)
under the European Unions Horizon 2020 research and innovation programme
(grant ERC CoG 725974). We thank Hideyuki Tachibana for the helpful
discussion. The contribution of Eliya Nachmani is part of a Ph.D. thesis
research conducted at Tel Aviv University.

## References

  * [1] Z.-Q. Wang, K. Tan, and D. Wang, “Deep learning based phase reconstruction for speaker separation: A trigonometric perspective,” in _ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2019, pp. 71–75.
  * [2] X. Li, X. Wu, and J. Chen, “A spectral-change-aware loss function for dnn-based speech separation,” in _ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2019, pp. 6870–6874.
  * [3] Z.-Q. Wang, J. L. Roux, D. Wang, and J. R. Hershey, “End-to-end speech separation with unfolded iterative phase reconstruction,” _arXiv preprint arXiv:1804.10204_ , 2018.
  * [4] Y. Luo and N. Mesgarani, “Tasnet: time-domain audio separation network for real-time, single-channel speech separation,” in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2018, pp. 696–700.
  * [5] E. Nachmani, Y. Adi, and L. Wolf, “Voice separation with an unknown number of multiple speakers,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 7164–7175.
  * [6] N. Zeghidour and D. Grangier, “Wavesplit: End-to-end speech separation by speaker clustering,” _arXiv preprint arXiv:2002.08933_ , 2020.
  * [7] D. Yu, M. Kolbæk, Z.-H. Tan, and J. Jensen, “Permutation invariant training of deep models for speaker-independent multi-talker speech separation,” in _2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2017, pp. 241–245.
  * [8] M. Kolbæk, D. Yu, Z.-H. Tan, and J. Jensen, “Multitalker speech separation with utterance-level permutation invariant training of deep recurrent neural networks,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 25, no. 10, pp. 1901–1913, 2017.
  * [9] N. Kanda, Y. Gaur, X. Wang, Z. Meng, and T. Yoshioka, “Serialized output training for end-to-end overlapped speech recognition,” _arXiv preprint arXiv:2003.12687_ , 2020.
  * [10] T. Ma, “Monaural source separation in the wild,” Ph.D. dissertation, Carnegie Mellon University Pittsburgh, PA, 2020.
  * [11] R. Martin and I. Cohen, “Single-channel speech presence probability estimation and noise tracking,” _Audio Source Separation and Speech Enhancement_ , pp. 87–106, 2018.
  * [12] O. Ernst, S. E. Chazan, S. Gannot, and J. Goldberger, “Speech dereverberation using fully convolutional networks,” in _2018 26th European Signal Processing Conference (EUSIPCO)_. IEEE, 2018, pp. 390–394.
  * [13] H. Erdogan, J. R. Hershey, S. Watanabe, and J. Le Roux, “Phase-sensitive and recognition-boosted speech separation using deep recurrent neural networks,” in _2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2015, pp. 708–712.
  * [14] E. Vincent, J. Barker, S. Watanabe, J. Le Roux, F. Nesta, and M. Matassoni, “The second ‘chime’speech separation and recognition challenge: Datasets, tasks and baselines,” in _2013 IEEE International Conference on Acoustics, Speech and Signal Processing_. IEEE, 2013, pp. 126–130.
  * [15] J. R. Hershey, Z. Chen, J. Le Roux, and S. Watanabe, “Deep clustering: Discriminative embeddings for segmentation and separation,” in _2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2016, pp. 31–35.
  * [16] Y. Isik, J. L. Roux, Z. Chen, S. Watanabe, and J. R. Hershey, “Single-channel multi-speaker separation using deep clustering,” _arXiv preprint arXiv:1607.02173_ , 2016.
  * [17] Z. Chen, Y. Luo, and N. Mesgarani, “Deep attractor network for single-microphone speaker separation,” in _2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2017, pp. 246–250.
  * [18] Y. Luo and N. Mesgarani, “Conv-tasnet: Surpassing ideal time–frequency magnitude masking for speech separation,” _IEEE/ACM transactions on audio, speech, and language processing_ , vol. 27, no. 8, pp. 1256–1266, 2019\.
  * [19] Y. Luo, Z. Chen, and T. Yoshioka, “Dual-path rnn: efficient long sequence modeling for time-domain single-channel speech separation,” _arXiv preprint arXiv:1910.06379_ , 2019.
  * [20] A. Shi, J. Han, and Z. Shi, “Toward the pre-cocktail party problem with tastas $+$,” _arXiv preprint arXiv:2009.03692_ , 2020.
  * [21] S. Markovich, S. Gannot, and I. Cohen, “Multichannel eigenspace beamforming in a reverberant noisy environment with multiple interfering speech signals,” _IEEE Transactions on Audio, Speech, and Language Processing_ , vol. 17, no. 6, pp. 1071–1086, 2009.
  * [22] T. N. Sainath, R. J. Weiss, K. W. Wilson, A. Narayanan, M. Bacchiani _et al._ , “Speaker location and microphone spacing invariant acoustic modeling from raw multichannel waveforms,” in _2015 IEEE Workshop on Automatic Speech Recognition and Understanding (ASRU)_. IEEE, 2015, pp. 30–36.
  * [23] Y. Luo, E. Ceolini, C. Han, S.-C. Liu, and N. Mesgarani, “Fasnet: Low-latency adaptive beamforming for multi-microphone audio processing,” _arXiv preprint arXiv:1909.13387_ , 2019.
  * [24] B. Laufer-Goldshtein, R. Talmon, and S. Gannot, “Global and local simplex representations for multichannel source separation,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 28, pp. 914–928, 2020.
  * [25] H. Tachibana, “Towards listening to 10 people simultaneously: An efficient permutation invariant training of audio source separation using sinkhorn’s algorithm,” in _2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021\.
  * [26] M. Yousefi, S. Khorram, and J. H. Hansen, “Probabilistic permutation invariant training for speech separation,” _arXiv preprint arXiv:1908.01768_ , 2019\.
  * [27] H. W. Kuhn, “The hungarian method for the assignment problem,” _Naval research logistics quarterly_ , vol. 2, no. 1-2, pp. 83–97, 1955.
  * [28] J. Munkres, “Algorithms for the assignment and transportation problems,” _Journal of the society for industrial and applied mathematics_ , vol. 5, no. 1, pp. 32–38, 1957.
  * [29] N. Tomizawa, “On some techniques useful for solution of transportation network problems,” _Networks_ , vol. 1, no. 2, pp. 173–194, 1971.
  * [30] J. Garofolo, D. Graff, D. Paul, and D. Pallett, “Csr-i (wsj0) complete ldc93s6a,” _Web Download. Philadelphia: Linguistic Data Consortium_ , vol. 83, 1993.
  * [31] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur, “Librispeech: an asr corpus based on public domain audio books,” in _2015 IEEE international conference on acoustics, speech and signal processing (ICASSP)_. IEEE, 2015, pp. 5206–5210.
  * [32] J. Cosentino, M. Pariente, S. Cornell, A. Deleforge, and E. Vincent, “Librimix: An open-source dataset for generalizable speech separation,” _arXiv preprint arXiv:2005.11262_ , 2020.
  * [33] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.


First Author et al.

mode = title]Electricity-gas integrated energy system optimal operation  in typical scenario of coal district considering hydrogen heavy trucks

1]Junjie Yin[orcid=0000-0001-7782-7274]

Conceptualization of this study, Methodology, Software,  Formal analysis, Investigation, Data Curation, Writing - Original Draft,Visualization
[1]School of Electrical Engineering, Southeast University, Nanjing 210096, China

1]Jianhua Wang[orcid=0000-0001-5185-4506]
[1]
Validation, Resources,  Writing - Review & Editing, Supervision, Project administration

1]Jun You
Writing - Review & Editing, Visualization, Supervision, Project administration

[cor1]Corresponding author

The coal industry contributes significantly to the social economy, but the emission of greenhouse gases puts huge pressure on the environment in the process of mining, transportation, and power generation. In the integrated energy system (IES), the current research about the power-to-gas (P2G) technology mainly focuses on the injection of hydrogen generated from renewable energy electrolyzed water into natural gas pipelines, which may cause hydrogen embrittlement of the pipeline and cannot be repaired. In this paper, sufficient hydrogen energy can be produced through P2G technology and coal-to-hydrogen (C2H) of coal gasification, considering the scenario of coal district is rich in coal and renewable energy. In order to transport the mined coal to the destination, hydrogen heavy trucks have a broad space for development, which can absorb hydrogen energy in time and avoid potentially dangerous hydrogen injection into pipelines and relatively expensive hydrogen storage. An optimized scheduling model of electric-gas IES is proposed based on second-order cone programming (SOCP). In the model proposed above, the closed industrial loop (including coal mining, hydrogen production, truck transportation of coal, and integrated energy systems) has been innovatively studied, to consume renewable energy and coordinate multi-energy. Finally, an electric-gas IES study case constructed by IEEE 30-node power system and Belgium 24-node natural gas network was used to analyze and verify the economy, low carbon, and effectiveness of the proposed mechanism.

* Advantages of hydrogen heavy trucks in typical scenario of coal industry district are discussed.

* The efficient closed loop of hydrogen energy from generation to utilization is proposed.

* P2G  including power-to-hydrogen and power-to-methane is developed, combined with coal gasification technology.

* The economy and low-carbon properties of the proposed IEGS  mechanism are analyzed.

Integrated energy systemPower-to-gas Coal-to-hydrogenHydrogen heavy trucksSecond-order cone programmingOptimized scheduling

§ INTRODUCTION

§.§ Motivations

As we all know, the coal district is a typical large energy load of industry. In the production process of mining, preparation, etc., coal mine machineries need to consume a lot of electricity and heat. In addition, heavy trucks transport coal to various destinations, which requires a large amount of energy to drive the engine and is essential in the coal industrial production process. Moreover, the coal industry areas usually suffer from  serious environmental pollution problems, the main reasons include but not limited to coal dust, factory emissions, vehicle exhaust, etc. Facing the tremendous environmental pressure and resource shortage, coal districts that rely on traditional energy supply methods are in urgent need of transformation.

Nowadays, the development of renewable energy, represented by wind power (WP)[A]WPWind power, solar photovoltaics (PV)[A]PVPhotovoltaics, etc., has become an inevitable trend. Generally speaking, coal resources are distributed in inland mountainous areas, and the local geographical features determine the abundant WP resources in these areas. However, the phenomenon of abandoning WP and PV is still serious, mainly due to the instability of renewable energy and the difficulty of being directly integrated into traditional energy networks.

Based on the above analysis, an energy mechanism for the typical coal districts has yet to be proposed. For coal districts, there are the following visions: giving full play to the advantages of rich wind energy in coal districts, designing practical and effective energy supply structure systems, and realizing smooth, low-carbon and environment-friendly transportation of coal resources. As an energy output terminal, it has important responsibilities for the entire energy network. Therefore, the following goals should be achieved: consuming renewable energy as much as possible, breaking the status quo of the separate design and independent operation of each existing energy supply system, and  establishing the integrated energy system (IES)[A]IESIntegrated energy system, in order to promote the optimization of the energy structure and finally achieve the grand goal of carbon neutrality.

§.§ Related Work

Regarding IES, many studies mainly focuse on the coupling and interaction of power systems and natural gas systems, , integrated electric-gas systems (IEGS)  [A]IEGSIntegrated electric-gas system[Zhang et al., 2018, Sheng et al., 2018]. According to the type of conversion gas, Power-to-gas (P2G) technology can be divided into two types: power-to-hydrogen (P2H) [A]P2HPower-to-hydrogenconversion and power-to-methane (P2M) [A]P2MPower-to-methaneconversion. The hydrogen and methane produced in P2G plants can be used in following ways [Gahleitner, 2013]:

* Further synthesis to methane or other hydrocarbon fuels through methanation equipment;

* Injection the mixture of hydrogen and methane into natural gas pipelines by hydrogen compressed natural gas (HCNG) technology [Zhou et al., 2020];

* Power generation with the help of internal combustion engine or combined heat and power (CHP) devices [Ni et al., 2018];

* Hydrogen storage into tanks after pressurization [Qiu et al., 2020];

* Hydrogen refueling stations for vehicles or the use of hydrogen in industry.

At present, some scholars are devoted to the research of hydrogen compressed natural gas (HCNG) equipment, and have achieved certain results in the study of the effect of hydrogen content percentage on integrated systems [Zhou et al., 2020]. However, due to the characteristics of low  density and high activity, blending hydrogen reduces the amount of energy delivered by the natural gas network under the same conditions. Hydrogen embrittlement may occur, which poses challenges to safety of  the overall system.  Due to the hazards of hydrogen embrittlement caused by the activity of hydrogen molecules, some studies are devoted to solving the location of hydrogen embrittlement cracks and performing subsequent repair work of steel pipelines through precise positioning [Chen et al., 2020]. Therefore, countries have established strict limits on hydrogen blending in natural gas networks (generally the volume of hydrogen is up to 6%) [IEA, 2020]. Then, the direct use of hydrogen converted by P2G technology in the CHP system involves multiple energy conversions, while each conversion is bound to bring about a large energy loss, so the gain is more than the loss. Also, High cost is the primary reason hindering large-scale storage of hydrogen energy, mainly due to pressurized systems and containers. Moreover, the current input cost of hydrogen buses in some cities is too high to make profits.

In a typical scenario of coal district, the conditions for the development of the hydrogen heavy trucks are unique: not only are hydrogen sources abundant, but also the application space of hydrogen heavy truck is broad. Different from the high hydrogen price in resource-deficient areas, the coal district has ample and cheap hydrogen, mainly from:

* More than 90% hydrogen is produced from fossil energy (, coal, natural gas) or alcohols, which emits greenhouse gases[Kunstle et al., 1978];

* Only 4% hydrogen is production obtained from renewable energy electrolysis, has more room for development[Pan et al., 2020]. (Relatively stable WP is more suitable than PV with greater daily volatility).

Due to the huge demand for coal transportation in these areas, the number of locally registered heavy trucks is also very large, usually reaching tens of thousands. The reasons and significance of developing hydrogen heavy trucks are as follows:

* Hydrogen heavy trucks can be directly developed, which are rigid demand, skipping low-profit fuel cell buses[Rolink and Rehtanz, 2013];

* Fossil energy consumption and greenhouse gas emissions can be reduced by replacing diesel heavy trucks;

* Compared with lithium power batteries, hydrogen heavy trucks are more suitable for heavy-duty and long-distance transportation, and have the advantages of longer cruising range, shorter charging time, and lighter weight;

* Solar and wind energy can be effectively absorbed, and hydrogen can be used as a carrier of energy capture, avoiding the negative impact of electricity generated by renewable energy directly connected to the grid [Wu et al., 2020];

* Electrolysis using renewable energy is usually the cleanest way, which can be considered as negative greenhouse gas production;

* The current hydrogen produced by electrolysis in various countries only accounts for 4-6% of the total hydrogen production, which is still the main improvement direction and policy preference area.

§.§ Contribution

[A]SOCPSecond-order cone programming
In order to solve the above-mentioned challenges, models of hydrogen heavy trucks and P2G equipment are established and  second-order cone programming (SOCP) are applied. The main contributions of this work are summarized below.

* The closed industrial loop between coal mining, hydrogen production, truck transportation of coal, and integrated energy systems has been innovatively proposed.

* P2G technology has been further expanded, by taking into account the detailed chemical reaction process of P2H and P2M. The proposed P2G technology adjusts electrolysis and methanation in real time to control the generation ratio of hydrogen and methane.

* Reasonable use of hydrogen energy, as the fuel of hydrogen heavy trucks, is proposed, which avoids hydrogen embrittlement and reduces the storage cost.

* In the scenario of coal district, the models of hydrogen heavy truck, P2G equipment and IEGS are innovatively proposed. Incorporating heavy trucks into the unified dispatch of IEGS is conducive to making full use of hydrogen energy, reducing carbon emissions, and achieving clean and environmental protection.

The remaining parts of this paper are organized as follows. The scenario of coal district considering hydrogen heavy trucks and P2G equipment model are introduced in Section II. Section III describes the problem formulation and transformation. Case study is provided to prove the validity of the proposed model in Section IV. Section V draws the conclusion on this paper.

§ TYPICAL SCENARIO OF COAL DISTRICT

Scenario of coal district considering hydrogen heavy trucks and P2G equipment model of WP  to H_2 and CH_4.

As shown in Fig. <ref>, the closed industrial loop between coal mining, hydrogen production, truck transportation of coal, and integrated energy systems has been innovatively proposed.

§.§ Mining, processing and transportation

§.§.§ Coal gasification

As summarized in section <ref>, more than 90% of hydrogen is produced by fossil energy, and the chemical reaction equation is as follows:

\begin{align}
&\ce{C + H2O ->[\triangle]CO  + H2}\ \Delta H=+131\  \rm{kJ/mol}\label{c1}\\
&\label{}\ce{CO + H2O ->CO2  + H2}\   \Delta H=-41\  \rm{kJ/mol}
\end{align}

Due to the scenario of
coal districts in this article, the process of producing hydrogen from fossil energy is only coal-to-hydrogen. Assuming that the weight of coal mined at the time slot $t$ is $M_t$ (ton), and $\beta_t$ (%)  is the ratio of mined coal used as the raw material for coal gasification, the amount of hydrogen produced by coal in this time slot is

\begin{equation}\label{}
f^{Coal}_{\ce{H2},t}=\alpha^{Coal}_t (\beta_t M_t)
\end{equation}

where $\alpha^{Coal}_t$ represents the efficiency of coal gasification at the  time slot $t$. This process is abbreviated as coal-to-hydrogen (C2H). Since the chemical reaction in (<ref>) requires heat absorption, the reaction equipment involved in the process is also one of the heat loads.
[A]C2HCoal-to-hydrogen of coal gasification
[V]$\alpha^{Coal}_t$Efficiency of coal gasification at the  time slot $t$
[V] $f^{Coal}_{\ce{H2},t}$Flow rate produced by C2H at the  time slot $t$
[V]$M_t$Weight of coal mined at the time slot $t$
[V]$\beta_t$Ratio of mined coal used as the raw material for coal gasification at the time slot $t$

§.§.§ Hydrogen heavy trucks

Comparison of heavy trucks with different power types.

Hydrogen truck    Electric vehicle(EV)   Diesel oil truck

Filling time    10-15 min    Several hours     10-15 min

Recharge mileage    500-750 mile    100-300 mile     500-750 mile

(Long-distance freight)   (Short-distance freight)

Impact on the grid    Buffer   Relatively large    No impact

Energy sustainability   Promising    Depend on battery technology     Large price fluctuations, limited reserves

Emission   Zero    Zero     High

Tab. <ref> shows the comparison between hydrogen heavy trucks and electric vehicles and diesel trucks in many aspects. It can be predicted that in the future when hydrogen energy production technology continues to increase and costs continue to fall, hydrogen heavy trucks will become more competitive.

In the actual production process, the amount of coal mined is directly proportional to the required transport capacity of heavy trucks. In this article, the load capacity and required quantity of each heavy truck are no longer considered separately. But from a global perspective, the relationship between the hydrogen consumption of hydrogen heavy trucks and the amount of coal mining is as follows:

\begin{equation}\label{}
f^{Truck}_{\ce{H2},t}=\alpha^{Truck}_t (1-\beta_t) M_t
\end{equation}

where $f^{Truck}_{\ce{H2},t}$ is the total hydrogen demand of heavy trucks at $t$ and $\alpha^{Truck}_t$ is the  transport hydrogen consumption coefficient of trucks. $(1-\beta_t) M_t$ represents the amount of coal that needs to be transported except for C2H.
[V] $f^{Truck}_{\ce{H2},t}$Total hydrogen demand of heavy trucks at $t$
[V]  $\alpha^{Truck}_t$Transport consumption coefficient of hydrogen trucks

§.§ P2G equipment model

§.§.§ Electrolysis & methanation technology

In the current field of industrial production of hydrogen, the main technology is proton exchange membrane (PEM) electrolysers. Electric conversion to hydrogen is achieved by electrolyzing water. The equation of the electrolyzed water reaction is shown in  Eq. (<ref>).

\begin{equation}\label{e18}
\ce{2H2O ->[{electrify}]2H2 ^ + O2 ^}
\end{equation}

Methanation is the conversion of COx to methane CH4 through hydrogenation. The reaction equations are shown in Eqs. (<ref>) & (<ref>) .

\begin{align}
&\label{e19}\ce{CO + 3H2 ->CH4  + H2O}\ \Delta H=-206\  \rm{kJ/mol}\\
&\label{e20}\ce{CO2 + 4H2 ->CH4  + 2H2O}\  \Delta H=-164\  \rm{kJ/mol}
\end{align}

§.§.§ P2H and P2M model

P2H and P2M technology have their own characteristics:

* P2H has higher reaction efficiency than P2M;

* P2H technology conversion conditions is less difficult than P2M , and the required cost is low;

* Hydrogen cannot be injected into existing natural gas pipelines on a large scale, while methane can;

* Methane is safer in view of the flammability limits[ISO 10156:2017(en), 2017] and explosion limits[ISO/IEC 80079-20-1:2017, 2017] ;

* The production and combustion of hydrogen do not involve carbon emissions, and CO2 is consumed during the production of methane;

* The calorific value of hydrogen is higher, which of methane is lower.

* Reaction efficiency constraints

\begin{align}
\begin{cases}
\eta_{elec,t}>\eta_{meth,t}\\
57\% \leq \eta_{elec} \leq 73\%\\
50\% \leq \eta_{meth} \leq 64\%
\end{cases}
\end{align}

[V]$\eta_{elec,t}$ Reaction efficiency of PEM electrolysis at $t$
[V]$\eta_{meth,t}$Reaction efficiency of methanation at $t$
where $\eta_{elec,t}$ and $\eta_{meth,t}$ are used to represent the reaction efficiency of the above two processes in the time slot $t$. Given two following assumptions:

(1) In each time slot, the conversion efficiency of the electrolytic cell and methanation equipment  $\eta_{elec,t}$, $\eta_{meth,t}$ remains unchanged, , fixed value;

(2) The chemical reaction process is a complete reaction and there is no reversible process.

* Conversion power consumption constraint

In the upper right corner of Fig. <ref>, the actual flow of the P2G equipment is drawn separately as Fig. <ref>, where each variable symbol (, flow rate and energy consumption parameters) is marked.
Schematic diagram of actual flow and variable symbols of P2G equipment.

\begin{align}
Cons_{\ce{H2},t}&=\alpha_{\ce{H2},t}  f_{\ce{H2},t} \\
Cons_{\ce{CH4},t}&=\alpha_{\ce{CH4elec},t} f^{'}_{\ce{H2},t}+\alpha_{\ce{CH4meth},t}  f_{\ce{CH4},t}\nonumber\\
&=4\alpha_{\ce{CH4elec},t} f_{\ce{CH4},t}+\alpha_{\ce{CH4meth},t}  f_{\ce{CH4},t}\label{e6}
\end{align}

[V]$Cons_{\ce{H2},t}$ Electricty consumption of producing H2
[V]$\alpha_{\ce{H2},t}$ Power consumption parameter corresponding to H2
[V]$f_{\ce{H2},t}$H2 flow rate produced by P2G equipment for hydrogen refueling station at  $t$
[V]$f^{'}_{\ce{H2},t}$H2 flow rate produced by P2G equipment for further methanation at  $t$

where $Cons_{\ce{H2},t}$  represents the electric energy consumption used by the P2G equipment to produce H2 for hydrogen refueling station; $f_{\ce{H2},t}$ is the H2 flow rate for hydrogen refueling station at the time slot $t$; $\alpha_{\ce{H2},t}$ is the power consumption parameter corresponding to H2 for hydrogen refueling station. Similarly, $Cons_{\ce{CH4},t}$ is the electric energy consumed to produce CH4;$f^{'}_{\ce{H2},t}$ is H2 flow rate produced for further methanation; $f_{\ce{CH4},t}$ is the flow rate of CH4; $\alpha_{\ce{CH4elec},t}$ represents the power consumption parameter  in process of electrolysis corresponding to CH4, and $\alpha_{\ce{CH4meth},t}$ is that in process of methanation. In order to simplify the situation, only the chemical reaction described in (<ref>) is considered here, that is, carbon dioxide and 4 times of the amount of the hydrogen produce methane under the action of high temperature and high pressure and a catalyst, , $f^{'}_{\ce{H2},t}=4 f_{\ce{CH4},t}$. Therefore the second equality is established in Eq. (<ref>).
[V]$\alpha_{\ce{CH4elec},t}$Power consumption parameter in process of electrolysis corresponding to CH4
[V]$\alpha_{\ce{CH4meth},t}$Power consumption parameter in process of methanation corresponding to CH4

* Flammability& explosion limit constraint[ISO 10156:2017(en), 2017, ISO/IEC 80079-20-1:2017, 2017]

The minimum concentration (%) of gas in the air under flammable conditions is defined as the Lower flammability limit (LFL), and the highest is the Upper flammability limit (UFL). Similarly, Explosion limits are expressed as lower explosive level (LEL) and Upper explosive limit (UEL), respectively.

\begin{align}
\rm{LFL}_{\ce{H2}} \leq \it{c}_{\ce{H2},t}\leq \rm{UFL}_{\ce{H2} }\\
\rm{LEL}_{\ce{H2}} \leq \it{c}_{\ce{H2},t} \leq \rm{UEL}_{\ce{H2} }
\end{cases}\\
\rm{LFL}_{\ce{CH4}} \leq \it{c}_{\ce{CH4},t} \leq \rm{UFL}_{\ce{CH4} }\\
\rm{LEL}_{\ce{CH4}} \leq \it{c}_{\ce{CH4},t} \leq \rm{UEL}_{\ce{CH4} }
\end{cases}
\end{align}

[V]$\it{c}_{\ce{H2}/\ce{CH4}}$  Concentration percentage of H2/CH4
[P]$\rm{LFL}$ Lower flammability limit of H2/CH4
[P]$\rm{UFL}$ Upper flammability limit of H2/CH4
[P]$\rm{LEL}$Lower explosive limit of H2/CH4
[P]$\rm{UEL}$Upper explosive limit of H2/CH4

* Calorific value constraint

The energy density of hydrogen is only about 25% of methane. When both release the same amount of energy, the volume of hydrogen is larger, and storage and transportation are more difficult. In this paper, $q_{\mathrm{H}_2}=119.96 $ MJ/kg, $q_{\mathrm{CH}_4}=50.00 $ MJ/kg.

§ PROBLEM FORMULATION

§.§ Minimize long-term total costs

\begin{align}
&\min _{\mathbf {P,F,\Pi}} C_{total}  =\sum_{i=1}^{N}\left(\sum_{i=1}^{T} C_{i}^{f}\left(P_{i, t}\right)+C_{i}^{U}+C_{i}^{D}\right)\\ \nonumber
& +\sum_{m=1}^{N_{g}} \sum_{t=1}^{T} C_{m, t}^{g}+\sum_{i=1}^{T} C_{{WP}, t}+\sum_{i=1}^{T} C_{{Truck}, t}-\sum_{i=1}^{T} E_{{Coal}, t}\\  \nonumber
&(3), (4), (8)-(12)  \nonumber\\
& \sum_{i=1}^{N_U} P_{i, t}=\sum_{i=1}^{N_{L}} P_{d, i, t}\label{power1}\\
& \sum P_{input,i,t}=\sum P_{output,i,t}\ \forall \textrm{power nodes}\\
& \sum_{i=1}^{N_U}\left(u_{i, t} P_{i, \max }-P_{i, t}\right) \geq \rho \sum_{i=1}^{N_{L}} P_{d, i, t} \\
& u_{i,t}=\left\{ \begin{array}{l}
\end{array} \right. \label{e11} \\
&u_{i, t} P_{i, \min } \leq P_{i} \leq u_{i, t} P_{i, \max } \\
& -R_{d} \leq P_{i,t}-P_{i, t-1} \leq R_{u}\label{e13}\\
&\sum_{k=t}^{t+\rm{T S}-1}(1-u_{i, t}) \geq \rm{T S}\it(u_{i, t\rm-1}-\it u_{i, t})\\
&\sum_{k=t}^{t+\rm{T O}-1} u_{i, k} \geq \rm{T O}\it\left(u_{i, t}-u_{i, t\rm-1}\right) \\
&C_{i, t}^{U} \geq \max\{H_{i}\left(u_{i, t}-u_{i, t-1}\right), 0\} \\
&C_{i, t}^{D} \geq \max\{J_{i}\left(u_{i, t-1}-u_{i, t}\right), 0\} \\
& P_{l, \min } \leq P_{l, t} \leq P_{l, \max }\label{power2}\\
& \sum_{n|\left( m,n \right) \in A}{f_{mn}}=\sum_{n|\left( n,m \right) \in A}{f_{nm}}+s_m\quad \forall m\in \mathcal{N}_g\label{heat1}\\
& \mathrm{sign}\left( f_{mn} \right) f_{mn}^{2}\le C_{mn}^{2}\left( p_{m}^{2}-p_{n}^{2} \right) \quad \forall m,n \in \mathcal{N}_g\label{e22}\\
& \mathrm{sign}\left( f_{mn} \right)=\left\{ \begin{array}{l}
1,\text{flow direction is positive}\\
\end{array} \right.\\
&\underline{s_m}\le s_m\le \overline{s_m}\quad \forall m\in \mathcal{N}_g \\
& \underline{p_m}\le p_m\le \overline{p_m}\quad \forall m\in \mathcal{N}_g\label{heat2}
\end{align}

[S]$\mathcal{N_U}$ Set of thermal units, $\mathcal{N_U}=\{1,2,\dots,i,\dots,N_U\}$
[S]$\mathcal{N_L}$ Set of load nodes, $\mathcal{N_L}=\{1,2,\dots,i,\dots,N_L\}$
[S]$\mathcal{N_g}$ Set of natural gas nodes, $\mathcal{N_g}=\{1,2,\dots,m,\dots,N_g\}$
[V]$P_{i, t}$  Output of thermal power unit $i$ at time $t$
[V]$P_{d, i, t}$ Load demand for electricity of node $j$ at time $t$
[V]$u_{i, t} $Start and stop status of unit $i$ at time $t$
[V]$P_{i, \max}$Maximum output value of unit $i$
[V]$P_{i, \min}$Minimum output value of unit $i$
[P]$\rho$Hot spare coefficient
[P]$R_{d}, R_{u}$ Unit down / up ramping speed
[P]TS, TO Minimum stop / start time
[P]$C_{i, t}^{U/D}$  Start/Stop cost of unit $i$ at time $t$
[P]$H_{i}, J_{i}$  Single start / stop cost of unit $i$
[V]$P_{l, t}$ Power flow of line $l$ at time $t$
[V]$f_{mn}$ Natural gas flow from node $m$ to node $n$
[V]$s_{m}$ Gas flow directly injected into the node $m$ from the source
[V]$p_m$Pressure value of node $m$
[P]$C_{mn}$ Weymouth constant
[V]$\underline{{s}_{m}},\overline{s_{m}}$ Threshold of natural gas source storage capacity
[V]$\underline{{p}_{m}},\overline{p_{m}}$Threshold of pipeline pressure
[V]$P_{input/output,i,t}$Power flow input/output of node $i$ at $t$

where the coal consumption cost function $C_{i}^{f}$ of thermal power unit $i$ can be expressed by the following quadratic equation:

\begin{align}
C_{i}^{f}\left( P_{i,t} \right) =a_iP_{i,t}^{2}+b_iP_{i,t}+c_i
\end{align}

where $C_{i}^{U}$ represents the start-up cost of unit $i$,
$C_{i}^{D}$ represents the shutdown cost of unit $i$,
$C_{m,t}^{g}$ represents the cost of the natural gas network.
For simplification, it is assumed that the punishment cost of abandoning WP ($C_{WP,t}$) is linearly positively related to the amount of electricity abandoning. $C_{{Truck}, t}$ indicates the cost of truck transportation, which is positively related to the load capacity of the truck, when the load does not exceed the load capacity of the truck. $E_{{Coal}, t}$ indicates the profit from selling coal, which is proportional to the weight of coal transportation.
[V]$C_{i}^{f}$Coal consumption cost of thermal power unit $i$
[V]$C_{m,t}^{g}$Cost of the natural gas network
[V]$C_{WP,t}$Punishment cost of abandoning WP
[V]$C_{{Truck}, t}$Cost of truck transportation
[V]$E_{{Coal}, t}$Profit from selling coal
Constraints Eqs. (<ref>)-(<ref>) are related to the power system, including power balance constraint, hot spare constraint, unit output constraint, unit ramp constraint, unit start and stop time constraint, start & stop cost constraint and line power flow safety constraint. And Eqs. (<ref>)-(<ref>) are related to the Natural gas system system, including flow balance constraint, Weymouth equations constraint[Borraz et al., 2016], source quantity constraint and pressure range constraint.

§.§ Problem transformation based on SOCP

In the unified power flow modeling and solving, the SOCP-based power flow model has been used in network planning  [Zhang et al., 2018] with binary decision variables, and the problems generated are all modeled as mixed integer SOCP (MISOCP) [Badesa et al., 2020].

In order to eliminate the non-linearity caused by pressure variables, make the following variable substitutions
\pi _m=p_{m}^{2}.%\quad \forall m\in \mathcal{N}g
[V]$\pi_{m}$ The square of the pressure value at node $m$
[V]$\underline{{\pi}_{m}},\overline{\pi_{m}}$Threshold of the square of the pressure value of $m$
Thus, the constraint (<ref>) can be rewritten as the following formula

\begin{equation}
\underline{\pi _m}\le \pi _m\le \overline{\pi _m}\quad \forall m\in \mathcal{N}g
\end{equation}

The non-convexity of the steady-state flow model is derived from Eq. (<ref>), where the absolute value
$\mathrm{sign}\left( f_{mn} \right) $ is non-smooth and non-differentiable. To solve this problem, a pair of binary variables $f_{mn}^+$ and $f_{mn}^-$  are introduced to represent the forward and backward flow directions of the pipe $m-n$, respectively. Therefore, Eq.<ref> is equivalently replaced with

\begin{align}\label{e31}
\left( \frac{ F_{m n} } {C_{m n}}\right)^{2}=\left(f_{m n}^{+}-f_{m n}^{-}\right)\left(\pi_{m}-\pi_{n}\right)
\forall m, n \in \mathcal{N}g
\end{align}

By using bilinear relaxation, the bilinear term on the right side of  (<ref>) can be replaced by  (<ref>)-(<ref>) (Note that since $f_{mn}^+$ and $f_{mn}^-$ are binary variables, this relaxation is strict):

\begin{align}\label{e32}
&\left( \frac{ F_{m n} } {C_{m n}}\right)^{2}=\lambda_{m n} \quad \forall m, n \in \mathcal{N}g\\\label{e33}
\lambda_{m n}  &\geq \pi_{n}-\pi_{m}+\left(f_{m n}^{+}-f_{m n}^{-}+1\right)\left(\pi_{m}^{l}-\pi_{n}^{u}\right) \\\label{e34}
\lambda_{m n} & \geq \pi_{n}-\pi_{m}+\left(f_{m n}^{+}-f_{m n}^{-}-1\right)\left(\pi_{m}^{u}-\pi_{n}^{l}\right) \\\label{e35}
\lambda_{m n}  &\geq \pi_{n}-\pi_{m}+\left(f_{m n}^{+}-f_{m n}^{-}+1\right)\left(\pi_{m}^{u}-\pi_{n}^{l}\right) \\\label{e36}
\lambda_{m n} & \geq \pi_{m}-\pi_{n}+\left(f_{m n}^{+}-f_{m n}^{-}-1\right)\left(\pi_{m}^{l}-\pi_{n}^{u}\right)
\end{align}

(<ref>) can be relaxed according to the cone format shown in (<ref>), where the standard SOC formula of (<ref>) is expressed as

\begin{align} \label{e37}
& \left(F_{m n} / C_{m n}\right)^{2} \leq \lambda_{m n} \quad \forall m, n \in \mathcal{N}g\\\label{e38}
&  \left\|\begin{array}{c}2 F_{m n} / C_{m n} \\ \lambda_{m n}-1\end{array}\right\|_{2} \leq \lambda_{m n}+1
\end{align}

§ CASE STUDY

IEGS case study diagram.

§.§ Software and solving tools

In order to verify the effectiveness of the method in this paper, the above optimization program was developed using MATLAB-YALMIP platform, and the CPLEX algorithm package was used to solve the nonlinear mixed integer programming problem. The hardware environment of the test system is Intel (R) Core (TM) i7-6500M CPU @ 2.50 GHz, 8 GB RAM, Win10 64 bit (operating system), Matlab R2019b (development environment), and the YALMIP version is R2020.

IBM CPLEX ILOG is a high-performance mathematical programming solver for linear programming, mixed integer programming.

The unit commitment problem is essentially a mixed integer programming (MIP) problem, which can be solved by MATLAB/CPLEX.

§.§ Initial data

Initial data of thermal power units.
2 *  Unit   2*Node     $P_{\max}$     $P_{\min}$   a   b   c   $R_u/R_d$   TS/TD   $H$   $J$

(p.u.)    (p.u.)     (ton/(p.u.)$^2$)    (ton/(p.u.))   (ton)   (p.u./h)   (h)   ($/time)   ($/time)

1    1    1.57     0.50   01524   38.5390   786.798   0.37   2   3937   19686

2    2    1.00    0.25    0.1058   46.1591   945.633   0.30   2   25000   12500

3    5    0.60     0.15   0.0280   40.3965   1049.998
0.15   2   15000   7500

4    8    0.80    0.20    0.0354   38.3055   1243.531   0.20   2   20000   10000

5    11    0.40    0.10   0.0211   36.327   1658.570   0.15   2   10000   5000

As shown in Fig. <ref>, the data of IEEE 30-node power system is utilized, and the rest of the data is shown in Tab. <ref>.  The IEEE 30-node test case represents a part of the American power system (located in the Midwest of the United States). The original IEEE 30-node scenario considers the processing conditions of 6 thermal power units. Here, Unit 6 connected to node 8 is changed to WP equipment.

The natural gas network data adopts Belgium 24-node network data [Borraz et al., 2016]. According to the internationally accepted natural gas sales rules, natural gas is priced according to its calorific value. The commonly used unit is Million Britain Thermal Unit (MBTU), while China's natural gas transaction is based on volume as a reference and the unit is cubic meter (m$^3$). It is known that 1 MBTU $\approx$ 28.3 m$^3$, and the unit conversion can be obtained: 1 $/MBTU $\approx$ 0.035 $/ m$^3$.

The maximum value of daily output power of the WP unit is set as 15,000 kW. The fluctuation data of WP output in each period is based on the project data of the University of Queensland in Australia. Now the unit output is increased by 30 times. Set the abandoning WP punishment parameter $\delta_{WP}$ = 0.08 $/kWh.
[P]$\delta_{WP}$ Abandoning WP punishment parameter

§.§ Simulation result

§.§.§ Typical daily operating cost

Typical daily WP output line chart.

Comparison of total operating costs with or without P2G.

Total operating costs     2 *  with P2G      2 * without P2G

$9.50763\times 10^6$   $10.30285\times 10^6$

The WP output data from February 27, 2019 to February 27, 2020 are selected, and the operation results in the proposed IEGS are shown in Fig. <ref>. The WP output at the same time basically conforms to winter>annual>summer, in which winter is selected from December to February of the following year, and summer is selected from June to August. Combined with the weather data measured by the wind project, WP output is influenced by natural factors such as temperature , wind speed and so on.

On the basis of ensuring the stable operation of WP output, in order to show the characteristics of energy-saving, clean, and fast absorption of P2G equipment, the comparison results obtained are shown in Tab. <ref>. Due to the independent operation of the power system and the natural gas system, and the direct connection of WP units to the grid, the electricity-gas system lacking a coupling and mutual aid relationship cannot absorb unstable and intermittent WP  in a timely manner. In particular, the cost of abandoning WP is penalized, leading to higher operating costs for the IEGS without P2G equipment. If the ecological benefits brought by the consumption of carbon dioxide gas by P2G equipment are taken into consideration, the advantages of P2G equipment are more prominent.

§.§.§ Low-carban

The energy conversion Sankey diagram of the proposed IEGS.

Select 14:00 data to draw Fig. <ref>, which is based on the annual WP output. It can be seen that among the initial energy sources, electricity and natural gas come from the power system and the natural gas system respectively, and renewable energy accounts for a considerable proportion. Through transitional energy conversion, the initial energy can be transformed into hydrogen (supplied to hydrogen heavy trucks), methane (can be directly injected into natural gas pipelines) and electricity (from power plants close to energy-consuming terminals).

In order to alleviate the greater environmental pressure and implement carbon neutral initiatives, carbon emissions trading is considered in this section[Ricke et al., 2018, Zhang et al., 2020]. According to the Status Report 2019 of International Carbon Action Partnership (ICAP), emissions trading systems of most countries in the world currently are in force, scheduled or under consideration. This article uses the carbon emissions trading of China and the European Union as a measurement standard, and calculates the operating costs of the three heavy trucks mentioned in subsection <ref>.

Comparison of total operating costs considering carbon emissions trading.

Total operating costs     2 * China      2 * European Union

Hydrogen truck   $9.47361\times 10^6$   $9.34901\times 10^6$

Electric vehicle   $9.56402\times 10^6$   $9.58834\times 10^6$

Diesel oil truck   $9.91750\times 10^6$   $10.20187\times 10^6$

According to the results of carbon emissions trading in 2020,  prices in China are about 30-40 ¥/tCO2, but that in the EU has exceeded 30 €/tCO2. Taking into account carbon emissions trading, the revised total cost is  ($C_{total}-C_{\ce{CO2}}$), where  $C_{\ce{CO2}}$ indicates the cost of CO2 emissions and proportional to the amount of CO2 emissions.
[V]$C_{\ce{CO2}}$Cost of CO2 emissions
Tab. <ref> shows that diesel oil trucks or LNG trucks will be levied high carbon emission fees, which will continue to increase as carbon neutrality advances. The current price of hydrogen heavy trucks is basically similar to that of electric vehicles, because in the proposed scenario, the methanation process can absorb a part of CO2. Therefore, hydrogen heavy trucks can not only achieve zero emissions on the basis of electric vehicles, but also absorb greenhouse gases, which are converted into profits here. With the continuous breakthrough of hydrogen energy technology and the inclination of government policies, the price of hydrogen energy will gradually drop and the scale of production will gradually expand. Further, hydrogen heavy trucks will have a broader prospect.

§.§.§ Unit output curve

Unit output curves under different hot spare coefficient.

When the hot spare coefficient $\rho$ is 0.05, the cumulative output ladder diagram of the unit is shown on the left side of Fig. <ref>-(a), which is basically in line with the distribution characteristics of power consumption peaks and valleys in reality, that is, there is a peak power consumption around noon during the day. In order to characterize the output status of each unit, draw the ladder diagram on the right side of Fig.<ref>-(a). Each curve is basically stable, and there are no several start and stop conditions within a short period of time, which plays a good role in the overall stability of the system.

Similarly, when the unit hot spare coefficient is 0.2, Fig. <ref>-(b) are available. The cumulative output ladder diagram of the unit is shown on the left side of Fig. <ref>-(b). The ladder diagram on the right side of Fig. <ref>-(b) shows the output status of each unit and the curves are basically stable. However, some problems can be found in comparison with $\rho$=0.05. The right side of Fig. <ref>-(b) shows that the output of one unit has a short start and stop near 19:00, which is caused by the higher hot spare. Under the premise that a certain unit operating life and system stability are lost, the system's ability to respond to emergencies can be greatly improved.

§.§.§ Sensitivity analysis of abandoning WP punishment parameters

In this section, adjust the abandoning WP punishment parameters appropriately, set to four values of 0.01, 0.05, 0.08, and 0.12 respectively. The operation results of IEGS are shown in Tab. <ref>. When the abandonment punishment parameter is too small, the operating cost is higher, second only to the operating cost under the condition of no P2G equipment, indicating that even if the punishment parameter is small, the amount of abandonment is large, which still leads to the system. The overall operating cost is high. With the increase of the abandonment punishment parameters, the operating cost of the system has been significantly reduced, which proves that considering the abandonment punishment cost in the optimization objective can make the system fully absorb PV power. However, if the punishment parameters for abandoning solar energy continue to increase, the operating cost of the system will increase. This is because a small amount of unconsumed PV will be multiplied by a larger punishment coefficient, which will cause a certain punishment cost.

Comparison of total operating costs with different abandoning PV punishment parameters.

Abandoning PV punishment parameters    Total operating costs

($/kWh)     ($/day )

0.01   $10.02773\times 10^6$

0.05   $9.54392\times 10^6$

0.08   $9.50763\times 10^6$

0.12   $9.84057\times 10^6$

Further, as shown in the dotted line in Fig. <ref>, it can be estimated that although there is a nonlinear relationship between the abandoning WP punishment parameter and the operating cost, there should be a minimum value as shown in Fig. <ref>. When the abandonment punishment parameter takes this value, the total operating cost of the system is the smallest, and the specific value needs to be further explored. Moreover, when the abandoning WP punishment parameter is less than the optimal value, the curve drops faster, and the difference between the values of the abandoning WP punishment parameter is not very large, but a large amount of abandonment energy multiplied by it will cause a larger punishment cost. Correspondingly, when the abandoning WP punishment parameter is greater than the optimal value, even if the abandoning WP punishment parameter becomes significantly larger, the cost curve rises slowly because the amount of abandonment remains at a small amount.

Curve line chart of the influence of $\delta_{WP}$ on $C_{total}$

§ CONCLUSION

In view of the gradual transformation of the energy structure and the continuous development of the Energy Internet, IESs are important research direction. This paper  draws the following conclusions:

* Typical daily operating cost is significantly reduced by 7.7%, by introducing the P2G model which converts seasonally volatile wind energy into hydrogen and methane.

* Hydrogen heavy trucks have certain advantages in the proposed coal districts scenario, , zero emissions and indirect consumption of renewable energy. Especially after the introduction of carbon emissions trading, the cost advantage is more prominent, compared with electric vehicles and diesel oil trucks.

* The proposed P2G equipment and C2H technology improve the ability to absorb WP, enhance the reliability of the system, and solve the problem of volatility. More importantly, the mechanism in the article circumvents the potential risks caused by hydrogen mixing in natural gas pipelines, the cost of which cannot be measured.

* The abandoning WP punishment parameter is considered in this article, and the optimal parameter is formulated according to the actual situation, which plays a certain guiding role for the future energy market leverage.

§ ACKNOWLEDGMENTS

This work was supported in part by the Natural Science Foundation of Jiangsu Province (BK20181283).

[Zhang et al., 2018]
authorY. Zhang, authorP. E. Campana,
authorY. Yang, authorB. Stridh,
authorA. Lundblad, authorJ. Yan,
titleEnergy flexibility from the consumer: Integrating
local electricity and heat supplies in a building,
journalApplied Energy volume223
(year2018) pages430–442.

[Sheng et al., 2018]
authorC. Sheng, authorZ. Wei, authorG. Sun,
authoret al.,
titleSteady state and transient simulation for
electricity-gas integrated energy systems by using convex optimisation,
journalIET Generation, Transmission & Distribution
volume12 (year2018) pages2199–2206.

[Gahleitner, 2013]
authorG. Gahleitner,
titleHydrogen from renewable electricity: An international
review of power-to-gas pilot plants for stationary applications,
journalInternational Journal of Hydrogen Energy
volume38 (year2013) pages2039–2061.

[Zhou et al., 2020]
authorS. Zhou, authorK. Sun, authorZ. Wu,
authoret al.,
titleOptimized operation method of small and medium-sized
integrated energy system for P2G equipment under strong uncertainty,
journalEnergy volume199
(year2020) pages117269.

[Ni et al., 2018]
authorL. Ni, authorW. Liu,
authorF. Wen, authorY. Xue,
authorZ. Dong, authorY. Zheng,
authorR. Zhang,
titleOptimal operation of electricity, natural gas and
heat systems considering integrated demand responses and diversified storage
journalJournal of Modern Power Systems and Clean Energy
volume6 (year2018) pages423–437.

[Qiu et al., 2020]
authorY. Qiu, authorS. Zhou, authorJ. Wang,
authorJ. Chou, authorY. Fang,
authorG. Pan, authorW. Gu,
titleFeasibility analysis of utilising underground
hydrogen storage facilities in integrated energy system: Case studies in
journalApplied Energy volume269
(year2020) pages115140.

[Chen et al., 2020]
authorY. Chen, authorH. Lu, authorJ. Liang,
authorA. Rosenthal, authorH. Liu,
authorG. Sneddon, authorI. McCarroll,
authorZ. Zhao, authorW. Li, authorA. Guo,
et al.,
titleObservation of hydrogen trapping at dislocations,
grain boundaries, and precipitates,
journalScience volume367
(year2020) pages171–175.

[IEA, 2020]
authorIEA, titleLimits on hydrogen blending in natural
gas networks, 2018, typeTechnical Report, IEA, Paris,

[Kunstle et al., 1978]
authorK. Kunstle, authorC. Koch,
authorK. Reiter, titleCoal gasification apparatus,
year1978. noteUnited States Patent No. 4095959,

[Pan et al., 2020]
authorG. Pan, authorW. Gu, authorH. Qiu,
authorY. Lu, authorS. Zhou, authorZ. Wu,
titleBi-level mixed-integer planning for
electricity-hydrogen integrated energy system considering levelized cost of
journalApplied Energy volume270
(year2020) pages115176.

[Rolink and Rehtanz, 2013]
authorJ. Rolink, authorC. Rehtanz,
titleLarge-scale modeling of grid-connected electric
journalIEEE Transactions on Power Delivery
volume28 (year2013) pages894–902.

[Wu et al., 2020]
authorX. Wu, authorH. Li,
authorX. Wang, authorW. Zhao,
titleCooperative operation for wind turbines and hydrogen
fueling stations with on-site hydrogen production,
journalIEEE Transactions on Sustainable Energy
volume11 (year2020) pages2775–2789.

[ISO 10156:2017(en), 2017]
authorISO 10156:2017(en), titleGas cylinders-Gases and
gas mixtures-Determination of fire potential and oxidizing ability for the
selection of cylinder valve outlets, typeTechnical Report, ISO,

[ISO/IEC 80079-20-1:2017, 2017]
authorISO/IEC 80079-20-1:2017, titleExplosive
atmospheres - Part 20-1: Material characteristics for gas and vapour
classification - Test methods and data, typeTechnical Report,
ISO/IEC, year2017.

[Borraz et al., 2016]
authorC. Borraz, authorR. Bent,
authorS. Backhaus, authorH. Hijazi,
authorP. V. Hentenryck,
titleConvex relaxations for gas expansion planning,
journalINFORMS Journal on Computing volume28
(year2016) pages645–656.

[Badesa et al., 2020]
authorL. Badesa, authorF. Teng,
authorG. Strbac,
titlePricing inertia and frequency response with diverse
dynamics in a mixed-integer second-order cone programming formulation,
journalApplied Energy volume260
(year2020) pages114334.

[Ricke et al., 2018]
authorK. Ricke, authorL. Drouet,
authorK. Caldeira, authorM. Tavoni,
titleCountry-level social cost of carbon,
journalNature Climate Change volume8
(year2018) pages895–900.

[Zhang et al., 2020]
authorX. Zhang, authorA. Löschel,
authorJ. Lewis, authorD. Zhang,
authorJ. Yan,
titleEmissions trading systems for global low carbon
energy and economic transformation,
journalApplied Energy volume279
(year2020) pages115858.