Next Article in Journal
HOME: 3D Human–Object Mesh Topology-Enhanced Interaction Recognition in Images
Next Article in Special Issue
Multi-Objective Material Logistics Planning with Discrete Split Deliveries Using a Hybrid NSGA-II Algorithm
Previous Article in Journal
Parabolic Hessian Equations Outside a Cylinder
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Optimal Control of PC-PLC Virus-Mutation and Multi-Delay Propagation Model in Distribution Network CPS

1
School of Mechanical and Electrical Engineering, Guangzhou University, Guangzhou 510006, China
2
Guangzhou Industry & Information Technology Institute for Intelligent Robotic Equipment, Guangzhou University, Guangzhou 510006, China
*
Authors to whom correspondence should be addressed.
Mathematics 2022, 10(16), 2840; https://doi.org/10.3390/math10162840
Submission received: 20 July 2022 / Revised: 7 August 2022 / Accepted: 8 August 2022 / Published: 9 August 2022

Abstract

:
The intelligent manufacturing of power systems has led to many challenges. The cyber-physical system (CPS) was introduced to solve the problem of insufficient integration of equipment and systems. It brings advantages, but also risks. In the distribution network CPS, malicious attacks on the PC-PLC communication network can cause significant incidents and affect system safety. The paper discusses two challenges, of possible mutated virus attacks and multi-delay in the PC-PLC coupled network. We present for the first time a virus-mutation and multi-delay propagation model. Then, to effectively control the virus propagation in the network and minimize the cost, the paper proposes five control measures, introduces their possible control combinations, and solves the optimal control problem with the Pontryagin maximum theorem. Finally, simulations verify the optimal control strategies for all combinations. By comparing the effects of maximum control, minimum control, average control, and optimal control, the optimal control strategy has been proven to be effective.

1. Introduction

With the development of intelligent manufacturing, the technology of electric power equipment has been improving, but there is still the problem of insufficient equipment and systems integration [1,2]. This problem leads to the intelligence of power production and operation not being fully realized in the smart grid [3]. The distribution network is one of the important links of power transmission, and intelligence is a major issue [4,5]. The cyber-physical system (CPS) [6] provides a new way to solve this problem. The information and physical space in the distribution network are organically integrated and deeply linked [7,8,9,10], forming the distribution network CPS. The effectiveness of resource allocation has increased with the addition of the cyber-physical system. Still, the emergence of new network dynamics characteristics has introduced different security risks to the distribution network CPS [11]. The secure operation of the distribution network faces a significant impact when the system is subject to active faults or malicious attacks [12]. In severe cases, it can even cause chain failures through information and physical interaction, threatening the safe operation of the distribution network and even the grid system. The distribution network CPS is usually based on a communication network between a PC and a PLC [13]. The protocols and devices used make them more vulnerable to cyber-attacks [14], leading to frequent security incidents. In January 2003, the Davis-Besse nuclear power plant in Ohio, USA, was attacked by the SQL Slammer worm, causing a system crash [15]. In June 2010, the malicious worm Stuxnet exploited vulnerabilities in industrial control systems to attack Iran’s nuclear power plants, destroying nearly one-fifth of the centrifuges [16]. In December 2015, the Ukrainian electricity network was struck by the BlackEnergy Trojan malware, resulting in a significant blackout of hundreds of thousands of households [17]. Therefore, distribution network system security has been particularly challenging [18].
Among the research methods for distribution network system security, Yang et al. [12] proposed a technique for CPS security assessment of distribution networks based on correlation matrix modeling analysis. Bai et al. [19] proposed a simulation method for distribution network faults based on CPS. Yao et al. [20] first proposed a mathematical model of PC-PLC worm reproduction based on an epidemiological kinetic approach. Among them, infectious epidemic dynamics play a guiding role in the transmission and treatment of diseases [21]. Since the transmission of computer viruses in particular networks is analogous to the spread of diseases, they are commonly exploited in industrial production [22,23]. As a result, by building a virus propagation model in computer networks, it is possible to examine the spread of viruses from their initial conditions, reveal propagation features, forecast change trends, and assess the effectiveness of various control methods [24,25]. Several studies have exploited infectious epidemic dynamics principles in the field of distribution network information security. Fu et al. [26] explored the dynamics of malware propagation by modeling epidemics with immunity and isolation. Sheng et al. [27] suggested an intelligent honeynet model to evaluate and prevent industrial virus proliferation in supervisory control and data acquisition (SCADA) networks.
Dynamical models are well known to be significantly impacted by delays. However, the above-proposed models ignore the delay in the virus propagation process. Zhang and Yang’s paper [28] mentioned the delay induced by the recovery node’s temporary immunity period. Wang and Chai [29] considered the delay caused by viral latency on top of it. Moreover, their study did not mention the isolation period of the quarantined nodes. In addition, the above delay exists during the virus propagation in the PC-PLC communication network. In this paper, we mainly consider the impact of the aforementioned delays on control.
Furthermore, there is a high likelihood of error codes being transmitted or the spread of the virus. This leads to the possibility of mutation of the virus [30]. Xu [31] mentioned the effect of mutated viruses during epidemic transmission in complex networks. Liu et al. [32] proposed a mutation model and an optimal control strategy for wireless rechargeable sensors (WRSN). In most networks, virus-mutations exist. Based on the above investigations, few people have studied the problem of optimal control under multi-delay and virus variation simultaneously in PC-PLC networks. Therefore, this paper focuses on this problem and provides optimal control strategies.
Using epidemic dynamics theory and optimal control theory to study the security of the distribution network CPS can reduce the harm of the virus invasion system with minimum cost. As a method to analyze optimal dynamic strategies [33], optimal control [34] is widely used in the safety assessment analysis of industrial systems. In the epidemic dynamics model, the primary means of virus control are vaccination, quarantine, and treatment [35]. The control methods in the PC-PLC coupled network described in this paper are mainly the installation of antivirus software or mitigating virus attacks by injecting patches [36].
In studying the virus propagation mechanism and optimal control strategies, this paper introduces a novel virus-mutation and multi-delay model to depict PC-PLC virus propagation. The contribution of this work is as follows:
  • Through introducing the effect of mutated viruses and delays into the existing PC-PLC model, a PLC-PC virus propagation model considering multi-delay is established.
  • Five model control measures are proposed, and 31 control combinations are given to solve the optimal control strategies by using the Pontryagin maximum theorem.
  • The optimal control outcomes of 31 control combinations of average constant and removing average constant control are described separately, and the optimal control strategies are studied and contrasted.
The subsequent sections of this paper are organized as follows: In Section 2, a PLC-PC virus-mutation and multi-delay propagation model is proposed. In Section 3, 31 control combinations of the model are proposed, and their optimal control strategies are solved. Section 4 shows the control simulation results, and the optimal control combination is analyzed. In Section 5, the conclusions of the proposed optimal control strategy are described.

2. Modeling Formulation

In this section, the problem of mutated viruses [37] is introduced into the PC-PLC coupled network, and the effect of multi-delay causes on the system is considered in a more realistic model for analysis.
For the problem of mutated viruses, parameters u A and u B are given as the probability of virus-mutation in PC and PLC, respectively. Mutated viruses are usually variable and unpredictable [38]. Therefore, mutated viruses should have higher infection rates, lower recovery rates, and isolation rates than viruses.
Delay has a significant impact on the stability of the system. Inspired by the work on multi-delay systems in [39,40], three different delay parameters, τ 1 , τ 2 , and τ 3 , are introduced in the model. Their definitions are as follows:
τ 1 : The delay between PC-exposed nodes and PLC-susceptible nodes transforming into infected nodes after being infected by PC is caused by virus latency time;
τ 2 : The delay of the quarantined nodes’ failure is caused by the isolation validity period;
τ 3 : The delay of the recovery nodes’ failure is caused by the immunity validity period.
In summary, this paper proposes a virus-mutation and multi-delay propagation model, and its differential equations are specified as follows:
{ d S A ( t ) d t = b A Λ A ( t ) + l A R A ( t τ 3 ) d A S A ( t ) β A S A ( t ) I 1 A ( t ) N A ( t ) d E A ( t ) d t = β A S A ( t ) I 1 A ( t ) N A ( t ) d A E A ( t ) ε E A ( t τ 1 ) d I 1 A ( t ) d t = ε E A ( t τ 1 ) + c 1 A Q A ( t τ 2 ) η 1 A I 1 A ( t ) d A I 1 A ( t ) u A I 1 A ( t ) γ 1 A I 1 A ( t ) d I 2 A ( t ) d t = u A I 1 A ( t ) + c 2 A Q A ( t τ 2 ) η 2 A I 2 A ( t ) d A I 2 A ( t ) γ 2 I 2 A ( t ) d Q A ( t ) d t = η 1 A I 1 A ( t ) + η 2 A I 2 A ( t ) c 1 A Q A ( t τ 2 ) c 2 A Q A ( t τ 2 ) d A Q A ( t ) ω A Q A ( t ) d R A ( t ) d t = γ 1 I 1 A ( t ) + γ 2 I 2 A ( t ) + ( 1 b A ) Λ A ( t ) + ω A Q A ( t ) l A R A ( t τ 3 ) d A R A ( t ) d S B ( t ) d t = Λ B ( t ) + l B R B ( t τ 3 ) + c B Q B ( t τ 2 ) η B S B ( t ) d B S B ( t ) β 1 B S B ( t ) I 1 B ( t ) N B ( t ) β 2 B S B ( t ) I 2 B ( t ) N B ( t ) β 1 A B S B ( t τ 1 ) I 1 A ( t τ 1 ) β 2 A B S B ( t τ 1 ) I 2 A ( t τ 1 ) d I 1 B ( t ) d t = β 1 B S B ( t ) I 1 B ( t ) N B ( t ) + β 1 A B S B ( t τ 1 ) I 1 A ( t τ 1 ) d B I 1 B ( t ) u B I 1 B ( t ) γ 1 B I 1 B ( t ) d I 2 B ( t ) d t = β 2 B S B ( t ) I 2 B ( t ) N B ( t ) + β 2 A B S B ( t τ 1 ) I 2 A ( t τ 1 ) u B I 1 B ( t ) d B I 2 B ( t ) γ 2 B I 2 B ( t ) d Q B ( t ) d t = η B S B ( t ) c B Q B ( t τ 2 ) d B Q B ( t ) ω B Q B ( t ) d R B ( t ) d t = ω B Q B ( t ) + γ 1 B I 1 B ( t ) + γ 2 B I 2 B ( t ) l B R B ( t τ 3 ) d B R B ( t ) ,
with initial conditions S A θ = 0 , I 1 A θ = 0 , I 2 A θ = 0 , R A θ = 0 , Q A θ = 0 , E A θ = 0 , S B θ = 0 , I 1 B θ = 0 , I 2 B θ = 0 , R B θ = 0 , Q B θ = 0 , θ max τ 1 , τ 2 , τ 3 ,   0 , S A 0 = Φ s A , I 1 A 0 = Φ I 1 A , S B 0 = Φ s B , I 2 A 0 = 0 , R A 0 = 0 , Q A 0 = 0 , E A 0 = 0 , I 1 B 0 = 0 , I 2 B 0 = 0 , R B 0 = 0 , Q B 0 = 0 , Φ s A , I 1 A 0 , Φ s B = Φ s A + Φ I 1 A .
Figure 1 shows the nodes’ state transformation process, and Table 1 lists its symbols and meanings. In Figure 1, the dashed line shows the delayed transformation.
For convenience and necessity, the model must satisfy the following conditions:
  • In the model, the total number of nodes at the PC and PLC ends is fixed (i.e., equipment inputs equal deaths). Thus, the PC network satisfies the constraint relation:
N A t = S A t + I 1 A t + I 2 A t + R A t + Q A t + E A t ,
the PLC network satisfies the constraint relation:
N B t = S B t + I 1 B t + I 2 B t + R B t + Q B t .
2.
Parameters ( β 1 A B , β 2 A B , β A , β 1 B , β 2 B , μ A , μ B , ε , b A , l A , l B , γ 1 A , γ 2 A , γ 1 B , γ 2 B , η 1 A , η 2 A , η B , c 1 A , c 2 A , c B , ω A , ω B , d A , d B ) in the model take values ranging from 0 to 1, except for Λ A , Λ B , τ 1 , τ 2   , and τ 3 .

3. Optimal Control

The optimal control strategy is established to make the distribution network cyber-physical system resistant to viruses and mutated viruses with minimal cost and loss.
To effectively combat the invasion of viruses and their mutated viruses, two measures are provided: injecting infected nodes and mutated infected nodes with immunity patches and injecting newly added PC nodes with pre-patch. Set γ 1 A = α 1 t , γ 2 A = α 2 t ,   γ 1 B = α 3 t ,   γ 2 B = α 4 t ,   1 b A = α 5 t . Among them, injecting immune patches of virus and mutated virus can improve the system’s recovery rate of infection and mutated infected nodes. Higher control intensity means a higher recovery rate of infected nodes and mutated infected nodes. Injecting newly added PC nodes with pre-patch aims to reduce the number of susceptible nodes in the system, thus reducing the rate of the virus and mutated virus propagation. Therefore, each control parameter is defined as follows:
α 1 t : The intensity of injected immune patch of virus in the PC nodes at the moment t ;
α 2 t : The intensity of injected immune patch of mutated virus in the PC nodes at the moment t ;
α 3 t : The intensity of injected immune patch of virus in the PLC nodes at the moment t ;
α 4 t : The intensity of injected immune patch of mutated virus in the PLC nodes at the moment t ;
α 5 t : The intensity of injected pre-patch in the newly added nodes at the moment t .
There are 31 combinations of the five control methods. This section discusses the optimal control strategies by comparing each combination’s control effectiveness.
First, the control strategies will be applied during 0 , t f . The control variable u t and the Lagrange function are introduced. The control variables u t for the optimal control strategy of each combination are given in Table 2.
Define the set of control function as:
U = α i t L 2 0 , t f : 0 t t f , 0 α i t 1 , i = 1 , 2 , 3 , 4 , 5 .
The objective function is given:
J ( u ) = 0 t f L ( u ( t ) ) d t = 0 t f E A t + I 1 A t + I 2 A t + I 1 B t + I 2 B t + m 2 α 1 2 t + n 2 α 2 2 t + p 2 α 3 2 t + q 2 α 4 2 t + s 2 α 5 2 t d t .
Our objective is to obtain the optimal control u *   such that J ( u * ) = m i n u ( u ) U J ( u ) . The optimal control strategy of each combination can be viewed as a particular case where some control parameters are set to constants. For example, when α 2 t = α 2 , α 3 t = α 3 ,   α 4 t = α 4 , α 5 t = α 5 , is the optimal control system expressed in Case 1. Parameters of optimal control for Cases 1–30 showed in Table 3.
The following is an example of Case 31. The optimal control objective function and the costing are as follows:
{ d S A t d t = 1 α 5 t Λ A t + l A R A t τ 3 d A S A t β A S A t I 1 A t N A t d E A t d t = β A S A ( t ) I 1 A ( t ) N A ( t ) d A E A ( t ) ε E A ( t τ 1 ) d I 1 A ( t ) d t = ε E A ( t τ 1 ) + c 1 A Q A ( t τ 2 ) η 1 A I 1 A ( t ) d A I 1 A ( t ) u A I 1 A ( t ) α 1 ( t ) I 1 A ( t ) d I 2 A ( t ) d t = u A I 1 A ( t ) + c 2 A Q A ( t τ 2 ) η 2 A I 2 A ( t ) d A I 2 A ( t ) α 2 ( t ) I 2 A ( t ) d Q A ( t ) d t = η 1 A I 1 A ( t ) + η 2 A I 2 A ( t ) c 1 A Q A ( t τ 2 ) c 2 A Q A ( t τ 2 ) d A Q A ( t ) ω A Q A ( t ) d R A ( t ) d t = α 1 ( t ) I 1 A ( t ) + α 2 ( t ) I 2 A ( t ) + α 5 ( t ) Λ A ( t ) + ω A Q A ( t ) l A R A ( t τ 3 ) d A R A ( t ) d S B ( t ) d t = Λ B ( t ) + l B R B ( t τ 3 ) + c B Q B ( t τ 2 ) η B S B ( t ) d B S B ( t ) β 1 B S B ( t ) I 1 B ( t ) N B ( t ) β 2 B S B ( t ) I 2 B ( t ) N B ( t ) β 1 A B S B ( t τ 1 ) I 1 A ( t τ 1 ) β 2 A B S B ( t τ 1 ) I 2 A ( t τ 1 ) d I 1 B ( t ) d t = β 1 B S B ( t ) I 1 B ( t ) N B ( t ) + β 1 A B S B ( t τ 1 ) I 1 A ( t τ 1 ) d B I 1 B ( t ) u B I 1 B ( t ) α 3 ( t ) I 1 B ( t ) d I 2 B ( t ) d t = β 2 B S B ( t ) I 2 B ( t ) N B ( t ) + β 2 A B S B ( t τ 1 ) I 2 A ( t τ 1 ) u B I 1 B ( t ) d B I 2 B ( t ) α 4 ( t ) I 2 B ( t ) d Q B ( t ) d t = η B S B ( t ) c B Q B ( t τ 2 ) d B Q B ( t ) ω B Q B ( t ) d R B ( t ) d t = ω B Q B ( t ) + α 3 ( t ) I 1 B ( t ) + α 4 ( t ) I 2 B ( t ) l B R B ( t τ 3 )
and
J α   1 , α   2 , α   3 , α   4 , α   5 = m i n u t U 0 t f E A t + I 1 A t + I 2 A t + I 1 B t + I 2 B t + m 2 α 1 2 t + n 2 α 2 2 t + p 2 α 3 2 t + q 2 α 4 2 t + s 2 α 5 2 t d t ,
where m , n , p , q , and s are the weight coefficients of the five control parameters, and m 2 α 1 2 t , n 2 α 2 2 t ,   p 2 α 3 2 t , q 2 α 4 2 t ,   s 2 α 5 2 t reflect the cost of the injected immune patch of virus in the PC node, the cost of the injected immune patch of mutated virus in the PC node, the cost of the injected immune patch of virus in the PLC node, the cost of the injected immune patch of mutated virus in the PLC node, and the cost of the injected pre-patch in the newly added nodes, respectively. Generally, they have the following limitations: 0 α i t 1 , i = 1 , 2 , 3 , 4 , 5 .
The goal of optimal control is to keep control costs as low as possible while reducing the number of infected and mutated infected nodes. Therefore, it is necessary to determine their optimal ratios under multiple control parameters. First, define the corresponding Hamiltonian function as:
H S A , E A , I 1 A , I 2 A , Q A , R A , S B , I 1 B , I 2 B , Q B , R B , α 1 , α 2 , α 3 , α 4 , α 5 , t = E A t + I 1 A t + I 2 A t + I 1 B t + I 2 B t + m 2 α 1 2 t + n 2 α 2 2 t + p 2 α 3 2 t + q 2 α 4 2 t + s 2 α 5 2 t + λ 1 t 1 α 5 t Λ A t + l A R A t τ 3 d A S A t β A S A t I 1 A t N A t + λ 2 t [ β A S A t I 1 A t N A t d A E A t ε E A t τ 1 ] + λ 3 t ε E A t τ 1 + c 1 A Q A t τ 2 η 1 A I 1 A t d A I 1 A t u A I 1 A t α 1 t I 1 A t + λ 4 t u A I 1 A t + c 2 A Q A t τ 2 η 2 A I 2 A t d A I 2 A t α 2 t I 2 A t + λ 5 t η 1 A I 1 A t + η 2 A I 2 A t c 1 A Q A t τ 2 c 2 A Q A t τ 2 d A Q A t ω A Q A t + λ 6 t α 1 t I 1 A t + α 2 t I 2 A t + α 5 t Λ A t + ω A Q A t l A R A t τ 3 d A R A t   + λ 7 t [ Λ B t + l B R B t τ 3 + c B Q B t τ 2 η B S B t d B S B t β 1 B S B t I 1 B t N B t β 2 B S B t I 2 B t N B t β 1 A B S B t τ 1 I 1 A t τ 1 β 2 A B S B t τ 1 I 2 A t τ 1 ] + λ 8 t β 1 B S B t I 1 B t N B t + β 1 A B S B t τ 1 I 1 A t τ 1 d B I 1 B t u B I 1 B t α 3 t I 1 B t + λ 9 t β 2 B S B t I 2 B t N B t + β 2 A B S B t τ 1 I 2 A t τ 1 + u B I 1 B t d B I 2 B t α 4 t I 2 B t   + λ 10 t η B S B t c B Q B t τ 2 d B Q B t ω B Q B t   + λ 11 t ω B Q B t + α 3 t I 1 B t + α 4 t I 2 B t l B R B t τ 3 d B R B t  
where   λ i t i = 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 , 11 is the covariate of the optimal control system.
According to the Pontryagin maximum theorem, the transversality conditions are:
λ i t f = 0 , i = 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 , 11 .
In addition to,
χ 0 , T τ i = 0 ,               t t f τ i , t f 1 ,                   t 0 , t f τ i , i = 1 , 2 , 3 .
The differential equations of covariates are as follows:
d λ 1 t d t = H S A t χ 0 ,   T τ 1 H S A τ 1 t + τ 1 χ 0 ,   T τ 2 H S A τ 2 t + τ 2 χ 0 ,   T τ 3 H S A τ 3 t + τ 3   = λ 1 t α 5 t 1 d A + d A + β A I 1 A N A I 1 A S A N A 2 λ 2 t β A I 1 A N A I 1 A S A N A 2 λ 6 t α 5 t d A ; d λ 2 t d t = H E A t χ 0 ,   T τ 1 H E A τ 1 t + τ 1 χ 0 ,   T τ 2 H E A τ 2 t + τ 2 χ 0 ,   T τ 3 H E A τ 3 t + τ 3   = 1 + λ 1 t α 5 t 1 d A + λ 2 t d A X 0 , T τ 1 ε λ 3 t + τ 1 λ 2 t + τ 1 λ 6 t α 5 t d A ; d λ 3 t d t = H I 1 A t χ 0 ,   T τ 1 H I 1 A τ 1 t + τ 1 χ 0 ,   T τ 2 H I 1 A τ 2 t + τ 2 χ 0 ,   T τ 3 H I 1 A τ 3 t + τ 3   = 1 + λ 1 t α 5 t 1 d A + β A S A N A I 1 A S A N A 2 λ 2 t β A S A N A I 1 A S A N A 2   + λ 3 t η 1 A + d A + u A + α 1 t λ 4 t u A λ 5 t η 1 A λ 6 t α 1 t + α 5 t d A   X 0 , T τ 1 λ 8 t + τ 1 β 1 A B S B λ 7 t + τ 1 β 1 A B S B ; d λ 4 t d t = H I 2 A t χ 0 ,   T τ 1 H I 2 A τ 1 t + τ 1 χ 0 ,   T τ 2 H I 2 A τ 2 t + τ 2 χ 0 ,   T τ 3 H I 2 A τ 3 t + τ 3   = 1 + λ 1 t α 5 t 1 d A + λ 4 t η 2 A + d A + α 2 t λ 5 t η 2 A λ 6 t α 2 t + α 5 t d A   + χ 0 , T τ 1 λ 7 t + τ 1 β 2 A B S B λ 9 t + τ 1 β 2 A B S B ;   d λ 5 t d t = H Q A t χ 0 ,   T τ 1 H Q A τ 1 t + τ 1 χ 0 ,   T τ 2 H Q A τ 2 t + τ 2 χ 0 ,   T τ 3 H Q A τ 3 t + τ 3   = λ 1 t α 5 t 1 d A + X 0 , T τ 2 c 1 A λ 5 t + τ 2 λ 3 t + τ 2 + X 0 , T τ 2 c 2 A λ 5 t + τ 2 λ 4 t + τ 2   + λ 5 t d A + ω A λ 6 t α 5 t d A + ω A ;   d λ 6 t d t = H R A t χ 0 ,   T τ 1 H R A τ 1 t + τ 1 χ 0 ,   T τ 2 H R A τ 2 t + τ 2 χ 0 ,   T τ 3 H R A τ 3 t + τ 3   = λ 1 t α 5 t 1 d A λ 6 t d A α 5 t 1 + X 0 , T τ 3 l A λ 6 t + τ 3 λ 1 t + τ 3 ;   d λ 7 t d t = H S B t χ 0 ,   T τ 1 H S B τ 1 t + τ 1 χ 0 ,   T τ 2 H S B τ 2 t + τ 2 χ 0 ,   T τ 3 H S B τ 3 t + τ 3   = λ 7 t η B + β 1 B I 1 B N B S B I 1 B N B 2 + β 2 B I 2 B N B S B I 2 B N B 2 + X 0 , T τ 1 λ 7 t + τ 1 β 1 A B I 1 A λ 8 t + τ 1 β 1 A B I 1 A   λ 8 t β 1 B I 1 B N B S B I 1 B N B + X 0 , T τ 1 λ 7 t + τ 1 β 2 A B I 2 A λ 9 t + τ 1 β 2 A B I 2 A   λ 9 t β 2 B I 2 B N B S B I 2 B N B 2 λ 10 t η B ;   d λ 8 t d t = H I 1 B t χ 0 ,   T τ 1 H I 1 B τ 1 t + τ 1 χ 0 ,   T τ 2 H I 1 B τ 2 t + τ 2 χ 0 ,   T τ 3 H I 1 B τ 3 t + τ 3   = 1 λ 7 t d B β 1 B S B N B I 1 B S B N B 2 λ 8 t β 1 B S B N B I 1 B S B N B 2 d B u B α 3 t λ 9 t u B λ 11 t α 3 t ;   d λ 9 t d t = H I 2 B t χ 0 ,   T τ 1 H I 2 B τ 1 t + τ 1 χ 0 ,   T τ 2 H I 2 B τ 2 t + τ 2 χ 0 ,   T τ 3 H I 2 B τ 3 t + τ 3   = 1 λ 7 t d B β 2 B S B N B I 2 B S B N B 2 λ 9 t β 2 B S B N B I 2 B S B N B 2 d B α 4 t λ 11 t α 4 t ;   d λ 10 t d t = H Q B t χ 0 ,   T τ 1 H Q B τ 1 t + τ 1 χ 0 ,   T τ 2 H Q B τ 2 t + τ 2 χ 0 ,   T τ 3 H Q B τ 3 t + τ 3   = λ 7 t d B + X 0 , T τ 2 c B λ 10 t + τ 2 λ 7 t + τ 2 + λ 10 t d B + ω B λ 11 t ω B ;   d λ 11 t d t = H R B t χ 0 ,   T τ 1 H R B τ 1 t + τ 1 χ 0 ,   T τ 2 H R B τ 2 t + τ 2 χ 0 ,   T τ 3 H R B τ 3 t + τ 3   = λ 7 t d B + X 0 , T τ 3 l B λ 11 t + τ 3 λ 7 t + τ 3 + λ 11 t d B .
Substituting the above conditions, we obtain:
H t α 1 α 1 t = α 1 * t = m α 1 * t + λ 3 t λ 6 t I 1 A * t ;   H t α 2 α 2 t = α 2 * t = n α 2 * t + λ 4 t λ 6 t I 2 A * t ;   H t α 3 α 3 t = α 3 * t = p α 3 * t + λ 8 t λ 11 t I 1 B * t ; H t α 4 α 4 t = α 4 * t = q α 4 * t + [ λ 9 t λ 11 t ] I 2 B * t ; H t α 5 α 5 t = α 5 * t = s α 5 * t + λ 1 t λ 6 t d A N A * t .
The optimal control ratio of the system is:
α 1 * t = m i n m a x 0 , λ 3 t λ 6 t I 1 A * t m , 1 ;   α 2 * t = m i n m a x 0 , λ 4 t λ 6 t I 2 A * t n , 1 ; α 3 * t = m i n m a x 0 , λ 8 t λ 11 t I 1 B * t p , 1 ; α 4 * t = m i n m a x 0 , [ λ 9 t λ 11 t ] I 2 B * t q , 1 ; α 5 * t = m i n m a x 0 , λ 1 t λ 6 t d A N A * t s , 1 .
Other combinations of optimal control solutions are unusual examples of the preceding solutions. Table 4 gives the specifics of these optimal solutions.

4. Numerical Simulations

In this section, we will simulate optimal control strategies for all combinations. Matlab is used in all of the simulations below. The parameters setting [21] is given in Table 5.
Patching a PLC requires exiting the node. Forced removal will affect access to information and network security. When a PC is removed, other nodes can usually take over [19]. Therefore, the PC will have less security impact and lower control costs. Mutated viruses are variable and unpredictable, causing high control costs. Thus, the value of cost weighting factors should follow: 0 < s < m < n < p < q . The cost weighting factors are set as follows: m = 300 ,   n = 400 ,   p = 500 ,   q = 600 ,   s = 1 . The simulation time is set as follows: t f = 50 . The delays are set as follows: τ 1 = 2 ,   τ 2 = 1.5 ,   τ 3 = 1 . The other initial conditions are set as follows: S A θ = 0 , I 1 A θ = 0 , I 2 A θ = 0 , R A θ = 0 , Q A θ = 0 , E A θ = 0 , S B θ = 0 , I 1 B θ = 0 , I 2 B θ = 0 , R B θ = 0 , Q B θ = 0 , θ 2 , 0 , S A 0 = 9000 , I 1 A 0 = 1000 , I 2 A 0 = 0 , R A 0 = 0 , Q A 0 = 0 , E A 0 = 0 , S B 0 = 10000 , I 1 B 0 = 0 , I 1 B 0 = 0 , R B 0 = 0 , Q B 0 = 0 .

4.1. Analysis of the Optimal Control Strategies with Average Constant Control

In Cases 1–30, the parameters are the average of the control parameters obtained in case 31, and are set as follows: α 1 = 0 t f   α 1 ( t ) / t f = 0 . 2217 , α 2 = 0 t f   α 2 ( t ) / t f = 0 . 0471 , α 3 = 0 t f   α 3 ( t ) / t f = 0 . 0445 , α 4 = 0 t f   α 4 ( t ) / t f = 0 . 0178 , α 5 = 0 t f   α 5 ( t ) / t f = 0 . 1370 .
The control parameters are shown in Table 6. The results of the optimal control strategies are shown in Figure 2, Figure 3, Figure 4 and Figure 5. As shown in Figure 2, Figure 3, Figure 4 and Figure 5, four small plots on the left show infected node curves over time (in clockwise order, the number of infected nodes in PC, mutated infected nodes in PC, infected nodes in PLC, and mutated infected nodes in PLC, respectively). In addition, the node states of the four nodes under average control, maximum control, minimum control, and optimal control are depicted. The last two graphs show optimal control parameter curves over time and control cost histograms for the four control scenarios. The control parameter tends to zero over time, as shown in the control parameter curve in Figure 2. The goal of this change is to lower the cost of control. It allows the system to have nearly the same virus suppression effect under optimal and maximum control at a significantly lower cost.
Table 6. Parameters of Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6.
Table 6. Parameters of Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6.
Case123456789101112131415
α 1 min00.22170.22170.22170.221700000.22170.22170.22170.22170.22170.2217
average0.22170.22170.22170.22170.22170.22170.22170.22170.22170.22170.22170.22170.22170.22170.2217
max10.22170.22170.22170.221711110.22170.22170.22170.22170.22170.2217
optimal-0.22170.22170.22170.2217----0.22170.22170.22170.22170.22170.2217
α 2 min0.047100.04710.04710.047100.04710.04710.04710000.04710.04710.0471
average0.04710.04710.04710.04710.04710.04710.04710.04710.04710.04710.04710.04710.04710.04710.0471
max0.047110.04710.04710.047110.04710.04710.04711110.04710.04710.0471
optimal0.0471-0.04710.04710.0471-0.04710.04710.0471---0.04710.04710.0471
α 3 min0.04450.044500.04450.04450.044500.04450.044500.04450.0445000.0445
average0.04450.04450.04450.04450.04450.04450.04450.04450.04450.04450.04450.04450.04450.04450.0445
max0.04450.044510.04450.04450.044510.04450.044510.04450.0445110.0445
optimal0.04450.0445-0.04450.04450.0445-0.04450.0445-0.04450.0445--0.0445
α 4 min0.01780.01780.017800.01780.01780.017800.01780.017800.017800.01780
average0.01780.01780.01780.01780.01780.01780.01780.01780.01780.01780.01780.01780.01780.01780.0178
max0.01780.01780.017810.01780.01780.017810.01780.017810.017810.01781
optimal0.01780.01780.0178-0.01780.01780.0178-0.01780.0178-0.0178-0.0178-
α 5 min0.13700.13700.13700.137000.13700.13700.137000.13700.137000.137000
average0.13700.13700.13700.13700.13700.13700.13700.13700.13700.13700.13700.13700.13700.13700.1370
max0.13700.13700.13700.137010.13700.13700.137010.13700.137010.137011
optimal0.13700.13700.13700.1370-0.13700.13700.1370-0.13700.1370-0.1370--
Case16171819202122232425262728293031
α 1 min0000000.22170.22170.22170.221700000.22170
average0.22170.22170.22170.22170.22170.22170.22170.22170.22170.22170.22170.22170.22170.22170.22170.2217
max1111110.22170.22170.22170.221711110.22171
optimal------0.22170.22170.22170.2217----0.2217-
α 2 min0000.04710.04710.04710000.04710000.047100
average0.04710.04710.04710.04710.04710.04710.04710.04710.04710.04710.04710.04710.04710.04710.04710.0471
max1110.04710.04710.04711110.04711110.047111
optimal---0.04710.04710.0471---0.0471---0.0471--
α 3 min00.04450.0445000.0445000.04450000.0445000
average0.04450.04450.04450.04450.04450.04450.04450.04450.04450.04450.04450.04450.04450.04450.04450.0445
max10.04450.0445110.0445110.04451110.0445111
optimal-0.04450.0445--0.0445--0.0445---0.0445---
α 4 min0.017800.017800.0178000.01780000.01780000
average0.01780.01780.01780.01780.01780.01780.01780.01780.01780.01780.01780.01780.01780.01780.01780.0178
max0.017810.017810.0178110.01781110.01781111
optimal0.0178-0.0178-0.0178--0.0178---0.0178----
α 5 min0.13700.137000.1370000.13700000.137000000
average0.13700.13700.13700.13700.13700.13700.13700.13700.13700.13700.13700.13700.13700.13700.13700.1370
max0.13700.137010.1370110.13701110.137011111
optimal0.13700.1370-0.1370--0.1370---0.1370-----
Case 31 contains the uncertainty of all control variables, and it should be optimal among all control combinations. For the rigor of the study, this still requires detailed simulation data to illustrate this conclusion. The implications of including different uncertainties in the cases are discussed below.
Cases 1–5 are the optimal control under one control measure. By observing Figure 3a–e, the impact of different control measures on infected system nodes is as follows:   α 1 t inhibits both infected nodes and mutated infected nodes of PC and PLC; α 2 t inhibits mutated infected nodes of PC and PLC;   α 3 t inhibits infected nodes and mutated infected nodes of PLC;   α 4 t inhibits mutated infected nodes of PLC;   α 5 t   slightly inhibits PC and PLC infected nodes and mutated infected nodes. The above is used as a basis to observe the control effect of the subsequent combinations.
First, we discuss the effect of the control of α 1 t and α 2 t (i.e., Figure 3f). Under this combination, the system effectively controls all infections. Under control with α 1 t and   α 3 t (i.e., Figure 3g), there is better control of the infection on the PLC side but less effective control of the mutated virus than the control of α 1 t and α 2 t . In control with α 1 t   and   α 4 t (i.e., Figure 3h), the system is only marginally more resistant to mutated viruses on the PLC side. When Cases 10–15 (i.e., Figure 3j and Figure 4a–e) are compared with Cases 22–25 and 30 (i.e., Figure 5b–e,j), it finds the system is difficult to suppress the spread of viruses and mutated viruses in the short term with the control without α 1 t .
According to Zhang’s paper [35], vaccination is not as effective as a treatment. It can significantly reduce the number of infections, but the number remains high in the short term. The purpose of vaccination is to reduce the susceptible population, similar to the parameter α 5 t . The control with α 5 t (i.e., Figure 3e) is effective in suppressing infected and mutated infected nodes but less effective in controlling the virus for a short period.
Following the above comparison, we found that effective virus control at the early stage can significantly reduce the rate of virus propagation. α 1 t is mainly used to control the virus that is put into the PC side system at the beginning. Therefore, α 1 t is the leading force in the five controls, and the rest are used as auxiliary means to cooperate with the controls and reduce costs. The effect of α 5 t is mainly in prevention, not removal. It effectively slows the spread of the virus, but relatively slowly. To show this effect, extending the system’s running time is necessary. Thus,   α 5 t   is a low-cost, long-term type of control.
Compared to combinations with   α 1 t (i.e., Figure 2, Figure 3f–i, Figure 4f–j, Figure 5a,f–i), there is no significant difference in effectiveness except cost. Thus, we must focus on their cost comparisons when considering these groups. Then, compared to typical combinations in Figure 5, the result still illustrates that α 1 t is critical and shows us that Case 31 is the least costly way to control. It also has the best control results by observing the magnified curve. It proves that Case 31 is the least costly and most effective control combination.

4.2. Analysis of the Optimal Control Strategies with Removing Average Constant Control

In this section, we remove the influence of the control constants (the averaged control amount is now taken as 0). The control parameters are shown in Table 7. This gives a more precise visualization of the impact of the optimal control parameters on the system. Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11 are the experiment’s results (i.e., Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6) after removing the control constants. Compared to the combination with α 1 t control (i.e., Figure 7, Figure 8f–i, Figure 9f–j, Figure 10a,f–i) and without α 1 t control (i.e., Figure 8b–e,j, Figure 9a–e, Figure 10b–e,j), other controls can significantly reduce the number of infected and mutated infected nodes when   α 1 t   control is not implemented. However, it is impossible to rapidly reduce the number of infected and mutated infected nodes to a low level in a short time.
In the comparison of Figure 6 and Figure 11, Cases 6, 26, and 31 are the more prominent. However, there is no significant difference in the effect, but Case 6 has a substantial difference in cost from the other two combinations. The remaining two groups are not remarkably different in cost and control effect, but Case 31 is slightly better than Case 26. In summary, among the 31 combinations of control strategies proposed, Case 31 is the optimal control combination.
To analyze the effectiveness of the optimal control strategy, the optimal control schemes are compared with the maximum, minimum, and average control parameters. In Figure 7, we can see that the effect of management is ordered as follows: maximum optimal > average > minimum , and the cost is ordered as follows: optimal < average < maximum < minimum . In maximum control, the control intensity is at its maximum all the time, so the effect of its control is slightly better than the effect of the optimal control. However, the cost is much higher than that of optimal control. Moreover, the average and minimum cases have poorer control effects and higher costs than the optimal. Therefore, it can be confirmed that the optimal control of Case 31 is effective.

5. Conclusions

To better control the propagation of viruses in the distribution network CPS, a novel model of virus propagation was developed after considering virus-mutation, immune delay, isolation delay, and infection delay. In addition, five control measures for the model are proposed. The 31 control strategies are listed, and the optimal control strategy under the optimal control combination is calculated.
In the simulation section, we compare the control effect and cost of all control combinations under maximum control, minimum control, average control, and optimal control. Then, the above situations are analyzed under average constant and without average constant control. It shows that injecting the immune patch of the virus into the PC network is the key to controlling the spread of viruses in the system. The optimal control under five control measures is derived as the optimal combination, and its effectiveness is verified.
This paper discusses a virus propagation model considering multi-delay and virus-mutation. Although optimal control is proposed, more research is needed to determine how to apply the proposed optimal control strategy to actual scenes. However, with the development of networks, control studies under mutation and multi-delay models are inevitable. We hope that the results of this paper can provide some reference for related researchers.

Author Contributions

Conceptualization, Y.S., Z.Q. and G.L.; methodology, Y.S., Z.Q. and G.L.; software, Y.S. and Z.Q.; validation, Y.S. and Z.Q.; formal analysis, Y.S. and Z.Q.; investigation, Y.S., G.L. and Z.L.; writing—original draft preparation, Y.S.; writing—review and editing, Y.S., G.L. and Z.L.; All authors have read and agreed to the published version of the manuscript.

Funding

The authors acknowledge funding received from the following science foundations: the National Natural Science Foundation of China (51975136, 52075109), the National Key Research and Development Program of China (2018YFB2000501), the Science and Technology Innovative Research Team Program in Higher Educational Universities of Guangdong Province (2017KCXTD025), the Special Research Projects in the Key Fields of Guangdong Higher Educational Universities (2019KZDZX1009), the Science and Technology Research Project of Guangdong Province (2017A010102014), the Industry-University-Research Cooperation Key Project of Guangzhou Higher Educational Universities (202235139), and the Guangzhou University Research Project (YJ2021002). They are all appreciated for supporting this work.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. He, B.; Bai, K.J. Digital twin-based sustainable intelligent manufacturing: A review. Adv. Manuf. 2021, 9, 1–21. [Google Scholar] [CrossRef]
  2. Jakobs, K. Standardisation of e-merging IoT-Applications: Past, Present and a Glimpse into the Future. In Proceedings of the 4th IEEE International Conference on Future Internet of Things and Cloud (FiCloud), Vienna, Austria, 22–24 August 2016. [Google Scholar]
  3. Alguliyev, R.; Imamverdiyev, Y.; Sukhostat, L. Cyber-physical systems and their security issues. Comput. Ind. 2018, 100, 212–223. [Google Scholar] [CrossRef]
  4. Zou, Y.; Chi, X.B.; Mai, H.B.; Zhang, S.Y. Discussion on Operation State Evaluation of Low Voltage Distribution Network. In Proceedings of the 4th International Conference on Intelligent Green Building and Smart Grid (IGBSG), Yichang, China, 6–9 September 2019; pp. 627–630. [Google Scholar]
  5. Haque, A.; Vo, T.H.; Nguyen, P.H. Distributed Intelligence: Unleashing Flexibilities for Congestion Management in Smart Distribution Networks. In Proceedings of the 4th IEEE International Conference on Sustainable Energy Technologies (ICSET), Hanoi, Vietnam, 14–16 November 2016; pp. 407–413. [Google Scholar]
  6. Wu, F.; Xu, L.L.; Li, X.; Kumari, S.; Karuppiah, M.; Obaidat, M.S. A Lightweight and Provably Secure Key Agreement System for a Smart Grid with Elliptic Curve Cryptography. IEEE Syst. J. 2019, 13, 2830–2838. [Google Scholar] [CrossRef]
  7. Zheng, B.W.; Yang, C.C.; Yang, J.; Bao, X.F. Analysis of Connection Modes Characteristics for Typical Distribution Networks with High Reliability. In Proceedings of the 3rd International Conference on Systems and Informatics (ICSAI), Shanghai, China, 19–21 November 2016; pp. 276–280. [Google Scholar]
  8. Chen, B.; Liu, Z.J.; Tang, Y.; Huang, J.Y.; Zhang, G.L.; Fan, Y.L. Typical Characteristics and Test Platform of CPS for Distribution Network. In Proceedings of the 7th IEEE Annual International Conference on CYBER Technology in Automation, Control, and Intelligent Systems (CYBER), Honolulu, HI, USA, 31 July–4 August 2017; pp. 1346–1350. [Google Scholar]
  9. Dong, W.J.; Liu, K.Y.; Hu, L.J.; Sheng, W.X.; Jia, D.L.; Deng, P. CPS Event Driving Method Based on Micro PMU of Distribution Network. In Proceedings of the 10th International Conference on Cyber-Enabled Distributed Computing and Knowledge Discovery (CyberC), Zhengzhou Univ, Zhengzhou, China, 18–20 October 2018; pp. 390–392. [Google Scholar]
  10. Zhao, Y.; Bai, M.K.; Liang, Y.; Ma, J.; Deng, P. Fault Modeling and Simulation Based on Cyber Physical System in Complex Distribution Network. In Proceedings of the China International Conference on Electricity Distribution (CICED), Tianjin, China, 17–19 September 2018; pp. 1566–1571. [Google Scholar]
  11. Zhou, X.; Yang, Z.; Ni, M.; Lin, H.S.; Li, M.L.; Tang, Y. Analysis of the Impact of Combined Information-Physical-Failure on Distribution Network CPS. IEEE Access 2020, 8, 44140–44152. [Google Scholar] [CrossRef]
  12. Yang, Z.; Zhou, X.; Ni, M.; Li, M.L.; Lin, H.S. A security Assessment of Distribution Network CPS Based on Association Matrix Modeling Analysis. In Proceedings of the 6th International Conference on Electrical Engineering, Control and Robotics (EECR)/3rd International Conference on Intelligent Control and Computing (ICICC), Xiamen, China, 10–12 January 2020. [Google Scholar]
  13. Osman, M.S.; Rahman, A.A.A.; Nordin, M.H. Development of User Interface for Cyber Physical System (CPS) Integrated With Material Handling System. In Proceedings of the Conference on Innovative Research and Industrial Dialogue (IRID), Melaka, Malaysia, 18 July 2018; pp. 170–171. [Google Scholar]
  14. Abe, S.; Tanaka, Y.; Uchida, Y.; Horata, S. Tracking Attack Sources based on Traceback Honeypot for ICS Network. In Proceedings of the 56th Annual Conference of the Society-of-Instrument-and-Control-Engineers-of-Japan (SICE), Kanazawa, Japan, 19–22 September 2017; pp. 717–723. [Google Scholar]
  15. Poulsen, K. Slammer Worm Crashed Ohio Nuke Plant Network. 2003. Available online: http://www.securityfocus.com/news/6767 (accessed on 5 July 2022).
  16. Farwell, J.P.; Rohozinski, R. Stuxnet and the Future of Cyber War. Survival 2011, 53, 23–40. [Google Scholar] [CrossRef]
  17. Cherepanov, A. BlackEnergy by the SSHBearDoor: Attacks against Ukrainian news media and electric industry. We Live Secur. 2016, 3, 2–4. [Google Scholar]
  18. Clark, A.; Li, Z.C.; Zhang, H.C. Control Barrier Functions for Safe CPS Under Sensor Faults and Attacks. In Proceedings of the 59th IEEE Conference on Decision and Control (CDC), Electr Network, Jeju, Korea, 14–18 December 2020; pp. 796–803. [Google Scholar]
  19. Bai, M.K.; Sheng, W.X.; Liang, Y.; Liu, K.Y.; Ye, X.S.; Kang, T.Y.; Wang, Y.K.; Zhao, Y. Research on Distribution Network Fault Simulation Based on Cyber Physics System. In Proceedings of the 3rd IEEE International Electrical and Energy Conference (CIEEC), Beijing, China, 7–9 September 2019; pp. 1112–1116. [Google Scholar]
  20. Yao, Y.; Sheng, C.; Fu, Q.; Liu, H.X.; Wang, D.J. A propagation model with defensive measures for PLC-PC worms in industrial networks. Appl. Math. Model. 2019, 69, 696–713. [Google Scholar] [CrossRef]
  21. Bjornstad, O.N.; Shea, K.; Krzywinski, M.; Altman, N. The SEIRS model for infectious disease dynamics. Nat. Methods 2020, 17, 557–558. [Google Scholar] [CrossRef]
  22. Lessler, J.; Azman, A.S.; Grabowski, M.K.; Salje, H.; Rodriguez-Barraquer, I. Trends in the Mechanistic and Dynamic Modeling of Infectious Diseases. Curr. Epidemiol. Rep. 2016, 3, 212–222. [Google Scholar] [CrossRef]
  23. Bailey, N.T. The Mathematical Theory of Infectious Diseases and Its Applications; Charles Griffin & Company Ltd.: High Wycombe Bucks, UK, 1975. [Google Scholar]
  24. Jia, Z.J.; Yang, Y.Y.; Guo, N. Research on Computer Virus Source Modeling with Immune Characteristics. In Proceedings of the 29th Chinese Control and Decision Conference (CCDC), Chongqing, China, 28–30 May 2017; pp. 4616–4619. [Google Scholar]
  25. Shah, D.; Zaman, T. Detecting Sources of Computer Viruses in Networks: Theory and Experiment. In Proceedings of the 2010 ACM SIGMETRICS International Conference on Measurement and Modeling of Computer Systems, New York, NY, USA, 14–18 June 2010; pp. 203–214. [Google Scholar]
  26. Fu, Q.; Yao, Y.; Sheng, C.; Yang, W. Interplay Between Malware Epidemics and Honeynet Potency in Industrial Control System Network. IEEE Access 2020, 8, 81582–81593. [Google Scholar] [CrossRef]
  27. Sheng, C.; Yao, Y.; Fu, Q.; Yang, W.; Liu, Y. Study on the intelligent honeynet model for containing the spread of industrial viruses. Comput. Secur. 2021, 111, 102460. [Google Scholar] [CrossRef]
  28. Zhang, Z.Z.; Yang, H.Z. Stability and Hopf Bifurcation in a Delayed SEIRS Worm Model in Computer Network. Math. Probl. Eng. 2013, 2013, 319174. [Google Scholar] [CrossRef]
  29. Wang, C.L.; Chai, S.X. Hopf bifurcation of an SEIRS epidemic model with delays and vertical transmission in the network. Adv. Differ. Equ. 2016, 2016, 100. [Google Scholar] [CrossRef]
  30. Avllazagaj, E.; Zhu, Z.Y.; Bilge, L.; Balzarotti, D.; Dumitras, T. When Malware Changed Its Mind An Empirical Study of Variable Program Behaviors in the RealWorld. In Proceedings of the 30th USENIX Security Symposium, Electr Network, Vancouver, BC, Canada, 11–13 August 2021; pp. 3487–3504. [Google Scholar]
  31. Xu, D.G.; Xu, X.Y.; Su, Z.F. Novel SIVR Epidemic Spreading Model with Virus Variation in Complex Networks. In Proceedings of the 27th Chinese Control and Decision Conference (CCDC), Qingdao, China, 23–25 May 2015; pp. 5164–5169. [Google Scholar]
  32. Liu, G.; Peng, Z.; Liang, Z.; Zhong, X.; Xia, X. Analysis and Control of Malware Mutation Model in Wireless Rechargeable Sensor Network with Charging Delay. Mathematics 2022, 10, 2376. [Google Scholar] [CrossRef]
  33. Zhu, L.H.; Wang, X.W.; Zhang, H.H.; Shen, S.L.; Li, Y.M.; Zhou, Y.D. Dynamics analysis and optimal control strategy for a SIRS epidemic model with two discrete time delays. Phys. Scr. 2020, 95, 035213. [Google Scholar] [CrossRef]
  34. Becerra, V.M. Solving optimal control problems with state constraints using nonlinear programming and simulation tools. IEEE Trans. Educ. 2004, 47, 377–384. [Google Scholar] [CrossRef]
  35. Yang, L.X.; Li, P.D.; Yang, X.F.; Xiang, Y.; Zhou, W.L. A Differential Game Approach to Patch Injection. IEEE Access 2018, 6, 58924–58938. [Google Scholar] [CrossRef]
  36. Xu, C.J.; Tang, X.H.; Liao, M.X.; He, X.F. Bifurcation analysis in a delayed Lokta-Volterra predator-prey model with two delays. Nonlinear Dyn. 2011, 66, 169–183. [Google Scholar] [CrossRef]
  37. Liu, G.Y.; Chen, J.Y.; Liang, Z.W.; Peng, Z.M.; Li, J.Q. Dynamical Analysis and Optimal Control for a SEIR Model Based on Virus Mutation in WSNs. Mathematics 2021, 9, 929. [Google Scholar] [CrossRef]
  38. Jia, X.Q.; Xiong, X.; Jing, J.W.; Liu, P. Using Purpose Capturing Signatures to Defeat Computer Virus Mutating. In Proceedings of the 6th International Conference on Information Security Practice and Experience, Seoul, Korea, 12–13 May 2010; pp. 153–171. [Google Scholar]
  39. Liu, G.Y.; Li, J.Q.; Liang, Z.W.; Peng, Z.M. Analysis of Time-Delay Epidemic Model in Rechargeable Wireless Sensor Networks. Mathematics 2021, 9, 978. [Google Scholar] [CrossRef]
  40. Zhang, L.; Liu, M.X.; Xie, B.L. Optimal control of an SIQRS epidemic model with three measures on networks. Nonlinear Dyn. 2021, 103, 2097–2107. [Google Scholar] [CrossRef]
Figure 1. The transformation process of nodes.
Figure 1. The transformation process of nodes.
Mathematics 10 02840 g001
Figure 2. Solution of the optimal control strategy for Case 31.
Figure 2. Solution of the optimal control strategy for Case 31.
Mathematics 10 02840 g002
Figure 3. Solutions of the optimal control strategies for Cases 1–10. (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4; (e) Case 5; (f) Case 6; (g) Case 7; (h) Case 8; (i) Case 9; (j) Case 10.
Figure 3. Solutions of the optimal control strategies for Cases 1–10. (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4; (e) Case 5; (f) Case 6; (g) Case 7; (h) Case 8; (i) Case 9; (j) Case 10.
Mathematics 10 02840 g003aMathematics 10 02840 g003b
Figure 4. Solutions of the optimal control strategies for Cases 11–20. (a) Case 11; (b) Case 12; (c) Case 13; (d) Case 14; (e) Case 15; (f) Case 16; (g) Case 17; (h) Case 18; (i) Case 19; (j) Case 20.
Figure 4. Solutions of the optimal control strategies for Cases 11–20. (a) Case 11; (b) Case 12; (c) Case 13; (d) Case 14; (e) Case 15; (f) Case 16; (g) Case 17; (h) Case 18; (i) Case 19; (j) Case 20.
Mathematics 10 02840 g004
Figure 5. Solutions of the optimal control strategies for Cases 21–30. (a) Case 21; (b) Case 22; (c) Case 23; (d) Case 24; (e) Case 25; (f) Case 26; (g) Case 27; (h) Case 28; (i) Case 29; (j) Case 30.
Figure 5. Solutions of the optimal control strategies for Cases 21–30. (a) Case 21; (b) Case 22; (c) Case 23; (d) Case 24; (e) Case 25; (f) Case 26; (g) Case 27; (h) Case 28; (i) Case 29; (j) Case 30.
Mathematics 10 02840 g005
Figure 6. Costs and effects of control combinations with better effects in Figure 2, Figure 3, Figure 4 and Figure 5.
Figure 6. Costs and effects of control combinations with better effects in Figure 2, Figure 3, Figure 4 and Figure 5.
Mathematics 10 02840 g006
Figure 7. Solution of the optimal control strategy for Case 31 (removing average constant).
Figure 7. Solution of the optimal control strategy for Case 31 (removing average constant).
Mathematics 10 02840 g007
Figure 8. Solutions of the optimal control strategies for Cases 1–10 (removing average constant). (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4; (e) Case 5; (f) Case 6; (g) Case 7; (h) Case 8; (i) Case 9; (j) Case 10.
Figure 8. Solutions of the optimal control strategies for Cases 1–10 (removing average constant). (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4; (e) Case 5; (f) Case 6; (g) Case 7; (h) Case 8; (i) Case 9; (j) Case 10.
Mathematics 10 02840 g008aMathematics 10 02840 g008b
Figure 9. Solutions of the optimal control strategies for Cases 11–20 (removing average constant). (a) Case 11; (b) Case 12; (c) Case 13; (d) Case 14; (e) Case 15; (f) Case 16; (g) Case 17; (h) Case 18; (i) Case 19; (j) Case 20.
Figure 9. Solutions of the optimal control strategies for Cases 11–20 (removing average constant). (a) Case 11; (b) Case 12; (c) Case 13; (d) Case 14; (e) Case 15; (f) Case 16; (g) Case 17; (h) Case 18; (i) Case 19; (j) Case 20.
Mathematics 10 02840 g009aMathematics 10 02840 g009b
Figure 10. Solutions of the optimal control strategies for Cases 21–30 (removing average constant). (a) Case 21; (b) Case 22; (c) Case 23; (d) Case 24; (e) Case 25; (f) Case 26; (g) Case 27; (h) Case 28; (i) Case 29; (j) Case 30.
Figure 10. Solutions of the optimal control strategies for Cases 21–30 (removing average constant). (a) Case 21; (b) Case 22; (c) Case 23; (d) Case 24; (e) Case 25; (f) Case 26; (g) Case 27; (h) Case 28; (i) Case 29; (j) Case 30.
Mathematics 10 02840 g010
Figure 11. Costs and effects of control combinations with better effects in Figure 7, Figure 8, Figure 9 and Figure 10.
Figure 11. Costs and effects of control combinations with better effects in Figure 7, Figure 8, Figure 9 and Figure 10.
Mathematics 10 02840 g011
Table 1. Definition of parameters in the system.
Table 1. Definition of parameters in the system.
NotationDefinition
S A t Number of PC susceptible nodes at time t
I 1 A t Number of PC infected nodes at time t
I 2 A t Number of PC mutated infected nodes at time t
R A t Number of PC recovered nodes at time t
Q A t Number of PC quarantined nodes at time t
E A t Number of PC exposed nodes at time t
S B t Number of PLC susceptible nodes at time t
I 1 B t Number of PLC infected nodes at time t
I 2 B t Number of PLC mutated infected nodes at time t
R B t Number of PLC recovered nodes at time t
Q B t Number of PLC quarantined nodes at time t
N A t Number of nodes in PC network at time t
N B t Number of nodes in PLC network at time t
Λ A t Number of PC nodes placements at time t
Λ B t Number of PLC nodes placements at time t
β 1 A B Infection rate of susceptible PLC nodes infected with the virus through PC nodes
β 2 A B Infection rate of susceptible PLC nodes infected with the mutated virus through PC nodes
β A Probability of exposure of susceptible PC nodes to the virus
β 1 B Infection rate of susceptible PLC nodes infected with the virus through PLC nodes
β 2 B . Infection rate of susceptible PLC nodes infected with the mutated virus through PLC nodes
μ A Mutaon rate of PC nodes infected with the virus
μ B Mutation rate of PLC nodes infected with the virus
ε Probability of transforming exposed nodes into infected nodes
b A Susceptibility rate of PC newly added nodes
l A Rate at which the recovered PC nodes lose immunity
l B Rate at which the recovered PLC nodes lose immunity
γ 1 A Recovered rate of PC infected nodes
γ 2 A Recovered rate of PC mutated infected nodes
γ 1 B Recovered rate of PLC infected nodes
γ 2 B Recovered rate of PLC mutated infected nodes
η 1 A Quarantined rate of PC infected nodes
η 2 A Quarantined rate of PC mutated infected nodes
η B Quarantined rate of infectious PC nodes
c 1 A Rate of infected PC nodes quarantined failure
c 2 A Rate of mutated infected PC nodes quarantined failure
c B Rate of susceptible PLC nodes quarantined failure
ω A Rate of quarantined PC nodes recovered
ω B Rate of quarantined PLC nodes recovered
d A Rate of PC nodes dead
d B Rate of PLC nodes dead
τ 1 Delay of Virus incubation period
τ 2 Delay of virus quarantined failure
τ 3 Delay of recovered nodes immune failure
Table 2. Control parameters, optimal control variables.
Table 2. Control parameters, optimal control variables.
CaseItems u t
1–5 α i , i = 1 , , n , n = 5 α i t , i = 1 , , n , n = 5
6–15 α i , α j , j = i + 1 , , n , i = 1 , , n 1 , n = 5 α i t , α j t , j = i + 1 , , n , i = 1 , , n 1 , n = 5
16–25 α i , α j , α k , k = j + 1 , , n , j = i + 1 , , n 1 , i = 1 , , n 2 , n = 5 α i t , α j t , α k t , k = j + 1 , , n , j = i + 1 , , n 1 , i = 1 , , n 2 , n = 5
26–30 α i , α j , α k , α l , l = k + 1 , , n , k = j + 1 , , n 1 , j = i + 1 , , n 2 , i = 1 , , n 3 , n = 5 α i t , α j t , α k t , α l t , l = k + 1 , , n , k = j + 1 , , n 1 , j = i + 1 , , n 2 , i = 1 , , n 3 , n = 5
31 α 1 , α 2 , α 3 , α 4 , α 5 α   1 t , α   2 t , α   3 t , α   4 t , α   5 t
Table 3. Optimal control parameters of Cases 1–30.
Table 3. Optimal control parameters of Cases 1–30.
Case123456789101112131415
α 1 α 1 t α 1 α 1 α 1 α 1 α 1 t α 1 t α 1 t α 1 t α 1 α 1 α 1 α 1 α 1 α 1
α 2 α 2 α 2 t α 2 α 2 α 2 α 2 t α 2 α 2 α 2 α 2 t α 2 t α 2 t α 2 α 2 α 2
α 3 α 3 α 3 α 3 t α 3 α 3 α 3 α 3 t α 3 α 3 α 3 t α 3 α 3 α 3 t α 3 t α 3
α 4 α 4 α 4 α 4 α 4 t α 4 α 4 α 4 α 4 t α 4 α 4 α 4 t α 4 α 4 t α 4 α 4 t
α 5 α 5 α 5 α 5 α 5 α 5 t α 5 α 5 α 5 α 5 t α 5 α 5 α 5 t α 5 α 5 t α 5 t
Case161718192021222324252627282930
α 1 α 1 t α 1 t α 1 t α 1 t α 1 t α 1 t α 1 α 1 α 1 α 1 α 1 t α 1 t α 1 t α 1 t α 1
α 2 α 2 t α 2 t α 2 t α 2 α 2 α 2 α 2 t α 2 t α 2 t α 2 α 2 t α 2 t α 2 t α 2 α 2 t
α 3 α 3 t α 3 α 3 α 3 t α 3 t α 3 α 3 t α 3 t α 3 α 3 t α 3 t α 3 t α 3 α 3 t α 3 t
α 4 α 4 α 4 t α 4 α 4 t α 4 α 4 t α 4 t α 4 α 4 t α 4 t α 4 t α 4 α 4 t α 4 t α 4 t
α 5 α 5 α 5 α 5 t α 5 α 5 t α 5 t α 5 α 5 t α 5 t α 5 t α 5 α 5 t α 5 t α 5 t α 5 t
Table 4. Solutions of optimal control problems 1–30.
Table 4. Solutions of optimal control problems 1–30.
C a s e 123456789101112131415
α 1 α 1 * t α 1 α 1 α 1 α 1 α 1 * t α 1 * t α 1 * t α 1 * t α 1 α 1 α 1 α 1 α 1 α 1
α 2 α 2 α 2 * t α 2 α 2 α 2 α 2 * t α 2 α 2 α 2 α 2 * t α 2 * t α 2 * t α 2 α 2 α 2
α 3 α 3 α 3 α 3 * t α 3 α 3 α 3 α 3 * t α 3 α 3 α 3 * t α 3 α 3 α 3 * t α 3 * t α 3
α 4 α 4 α 4 α 4 α 4 * t α 4 α 4 α 4 α 4 * t α 4 α 4 α 4 * t α 4 α 4 * t α 4 α 4 * t
α 5 α 5 α 5 α 5 α 5 α 5 * t α 5 α 5 α 5 α 5 * t α 5 α 5 α 5 * t α 5 α 5 * t α 5 * t
C a s e 161718192021222324252627282930
α 1 α 1 * t α 1 * t α 1 * t α 1 * t α 1 * t α 1 * t α 1 α 1 α 1 α 1 α 1 * t α 1 * t α 1 * t α 1 * t α 1
α 2 α 2 * t α 2 * t α 2 * t α 2 α 2 α 2 α 2 * t α 2 * t α 2 * t α 2 α 2 * t α 2 * t α 2 * t α 2 α 2 * t
α 3 α 3 * t α 3 α 3 α 3 * t α 3 * t α 3 α 3 * t α 3 * t α 3 α 3 * t α 3 * t α 3 * t α 3 α 3 * t α 3 * t
α 4 α 4 α 4 * t α 4 α 4 * t α 4 α 4 * t α 4 * t α 4 α 4 * t α 4 * t α 4 * t α 4 α 4 * t α 4 * t α 4 * t
α 5 α 5 α 5 α 5 * t α 5 α 5 * t α 5 * t α 5 α 5 * t α 5 * t α 5 * t α 5 α 5 * t α 5 * t α 5 * t α 5 * t
Table 5. Simulation parameters.
Table 5. Simulation parameters.
NotationDefinition
β 1 A B 0.0000008
β 2 A B 0.0000016
β A 0.08
β 1 B 0.000008
β 2 B 0.00001
μ A 0.005
μ B 0.001
ε 0.05
l A 0.001
l B 0.001
η 1 A 0.001
η 2 A 0.0008
c 1 A 0.001
c 2 A 0.004
c B 0.004
ω A 0.001
ω B 0.001
d A 0.004
d B 0.004
Table 7. Parameters of Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11.
Table 7. Parameters of Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11.
Case123456789101112131415
α 1 min000000000000000
average0.221700000.22170.22170.22170.2217000000
max100001111000000
optimal-0000----000000
α 2 min000000000000000
average00.04710000.04710000.04710.04710.0471000
max010001000111000
optimal0-000-000---000
α 3 min000000000000000
average000.04450000.0445000.0445000.04450.04450
max001000100100110
optimal00-000-0.04450-00--0
α 4 min000000000000000
average0000.01780000.0178000.017800.017800.0178
max000100010010101
optimal000-000-00-0-0-
α 5 min000000000000000
average00000.13700000.1370000.137000.13700.1370
max000010001001011
optimal0000-000-00-0--
Case16171819202122232425262728293031
α 1 min0000000000000000
average0.22170.22170.22170.22170.22170.221700000.22170.22170.22170.221700.2217
max1111110000111101
optimal------0000----0-
α 2 min0000000000000000
average0.04710.04710.04710000.04710.04710.047100.04710.04710.047100.04710.0471
max1110001110111011
optimal---000---0---0--
α 3 min0000000000000000
average0.0445000.04450.044500.04450.044500.04450.04450.044500.04450.04450.0445
max1001101101110111
optimal-00--0--0---0---
α 4 min0000000000000000
average00.017800.017800.01780.017800.01780.01780.017800.01780.01780.01780.0178
max0101011011101111
optimal0-0-0--0---0----
α 5 min0000000000000000
average000.137000.13700.137000.13700.13700.137000.13700.13700.13700.13700.1370
max0010110111011111
optimal00-0--0---0-----
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Su, Y.; Qiu, Z.; Liu, G.; Liang, Z. Optimal Control of PC-PLC Virus-Mutation and Multi-Delay Propagation Model in Distribution Network CPS. Mathematics 2022, 10, 2840. https://doi.org/10.3390/math10162840

AMA Style

Su Y, Qiu Z, Liu G, Liang Z. Optimal Control of PC-PLC Virus-Mutation and Multi-Delay Propagation Model in Distribution Network CPS. Mathematics. 2022; 10(16):2840. https://doi.org/10.3390/math10162840

Chicago/Turabian Style

Su, Yingying, Zijing Qiu, Guiyun Liu, and Zhongwei Liang. 2022. "Optimal Control of PC-PLC Virus-Mutation and Multi-Delay Propagation Model in Distribution Network CPS" Mathematics 10, no. 16: 2840. https://doi.org/10.3390/math10162840

APA Style

Su, Y., Qiu, Z., Liu, G., & Liang, Z. (2022). Optimal Control of PC-PLC Virus-Mutation and Multi-Delay Propagation Model in Distribution Network CPS. Mathematics, 10(16), 2840. https://doi.org/10.3390/math10162840

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop