SYSTEMS AND METHODS FOR ASSESSING THE VULNERABILITY OF CYBER-PHYSICAL SYSTEMS

Abstract

Methods for securing cyber-physical systems include attack generation systems, methods, and devices configured to assess vulnerabilities of cyber-physical systems include an example method of training an attack generative model using an attack policy to generate an attack dataset, training discriminators using a random attack dataset and generated attack dataset and training a generator using the trained discriminators.

Description

BACKGROUND

Industrial and infrastructure systems commonly include supervisory, control, and data acquisition (SCADA) systems. SCADA includes data collection, estimation, and/or decision processes for industrial facilities. SCADA systems are commonly implemented using networks that can include any number of computing devices, sensors, user controls, actuators, etc. SCADA networks can sometimes include internet access, and/or can cover large areas. As a result, SCADA systems can be difficult to secure, and can be significant vulnerabilities for industry and infrastructure installments. There are benefits to improving cybersecurity for such facilities, in particular the networks and devices used to monitor and/or operate infrastructure and industrial systems.

SUMMARY

In some aspects, implementations of the present disclosure include a computer-implemented method of training a machine learning model, the method including: generating a random parameter vector; generating a random attack dataset; inputting a plurality of samples of the random attack dataset into a generator; generating, by the generator, a generated attack dataset; training a discriminator using the random attack dataset and the generated attack dataset; and training the generator using the trained discriminators and a loss function.

In some aspects, implementations of the present disclosure include a computer-implemented method, wherein the generated attack dataset includes an attack policy.

In some aspects, implementations of the present disclosure include a computer-implemented method, wherein the attack policy includes a ramp attack.

In some aspects, implementations of the present disclosure include a computer-implemented method, wherein the attack policy includes a sensor attack.

In some aspects, implementations of the present disclosure include a computer-implemented method, wherein the generator includes a deep neural network.

In some aspects, implementations of the present disclosure include a computer-implemented method, wherein the generator is trained with the loss function.

In some aspects, implementations of the present disclosure include a computer-implemented method for training a machine learning model, the method including: receiving a generative model parameter vector; generating sample points from a prior distribution of the generative model parameter vector; generating attack policy parameters simulating simulated attack policy parameters by a system simulation training a generative model by the simulated attack policy parameters.

In some aspects, implementations of the present disclosure include a computer-implemented method, wherein training the generative model includes a deep neural network

In some aspects, implementations of the present disclosure include a computer-implemented method, wherein training the generative model includes selecting a best loss value by.

In some aspects, implementations of the present disclosure include a computer-implemented method, wherein the attack policy parameters include a ramp attack policy.

In some aspects, implementations of the present disclosure include a computer-implemented method, wherein the attack policy parameters include a sine attack policy.

In some aspects, implementations of the present disclosure include a computer-implemented method, wherein the attack policy parameters include a pulse attack policy.

In some aspects, implementations of the present disclosure include a cybersecurity system including: a physical plant a network configured for communications and/or control of the physical plant; a cybersecurity controller operably connected to the network, wherein the cybersecurity controller includes a processor and memory with instructions stored thereon, that, when executed by the processor cause the processor to: receive a trained attack generative model; simulate, by the trained attack generative model, a plurality of attacks on the physical plant; determine an attack effectiveness of each of the plurality of attacks; and control access to the network based on the effectiveness of each of the plurality of attacks.

In some aspects, implementations of the present disclosure include a system, wherein the physical plant includes a networked pipeline system.

In some aspects, implementations of the present disclosure include a system, wherein the networked pipeline system includes a plurality of pressure sensors operably coupled to the network, and controlling access to the network includes securing at least one of the plurality of pressure sensors.

In some aspects, implementations of the present disclosure include a system, wherein the physical plant includes a power grid.

In some aspects, implementations of the present disclosure include a system, wherein the power grid includes a plurality of meters configured to measure the power and frequency of the power grid, and controlling access to the network includes securing at least one of the plurality of meters.

In some aspects, implementations of the present disclosure include a system, wherein the physical plant includes a cyber-physical system.

In some aspects, implementations of the present disclosure include a system, wherein the physical plant includes an industrial facility.

In some aspects, implementations of the present disclosure include a system, wherein the controller further contains instructions that cause the processor to simulate a vulnerability of the physical plant and control the physical plant based on the vulnerability.

In some aspects, implementations of the present disclosure include a system, wherein the controller is further configured to control the network device based on the sampled actions implemented on the system model.

Other systems, methods, features and/or advantages will be or may become apparent to one with skill in the art upon examination of the following drawings and detailed description. It is intended that all such additional systems, methods, features and/or advantages be included within this description and be protected by the accompanying claims.

BRIEF DESCRIPTION OF THE DRAWINGS

The components in the drawings are not necessarily to scale relative to each other. Like reference numerals designate corresponding parts throughout the several views.

FIG. 1 illustrates a cybersecurity system, according to implementations of the present disclosure.

FIG. 2A illustrates a method of training a machine learning model to evaluate a vulnerability of a physical system, according to implementations of the present disclosure.

FIG. 2B illustrates a method of training a machine learning model to evaluate a vulnerability of a physical system, according to implementations of the present disclosure.

FIG. 3A illustrates an example implementation of the method of FIG. 2A.

FIG. 3B illustrates an example implementation of the method of FIG. 2B.

FIG. 4 illustrates an example computing device.

FIG. 5 illustrates a schematic diagram of SCADA for networked pipeline systems under adversarial attacks.

FIG. 6 illustrates a schematic description of an example implementation of a vulnerability analysis model, according to implementations of the present disclosure.

FIG. 7 illustrates a four-node pipeline network simulation platform used in a study of an example implementation of the present disclosure.

FIG. 8 illustrates an example of nominal control performance with disturbance injection between 100 s and 120 s.

FIG. 9 illustrates a study of an example implementation of an attack generator in a system simulation.

FIG. 10 illustrates loss curves of example discriminators at different training epochs in a study of an example implementation of the present disclosure.

FIG. 11 illustrates example control performance under generated attacks from a trained attack generator in a study of an example implementation of the present disclosure.

FIG. 12 illustrates an example BDD residual under generated attacks, in a study of an example implementation of the present disclosure

FIG. 13 illustrates an example control system of a networked cyber-physical system (NCPS) according to implementations of the present disclosure.

FIG. 14 illustrates an example training framework of an attack generative model with pre-defined attack policy and custom loss function, according to implementations of the present disclosure.

FIG. 15 illustrates example attack policies used in an example implementation of the present disclosure.

FIG. 16 illustrates a schematic description of a cyber-physical IEEE 14-bus system, according to an example implementation of the present disclosure.

FIG. 17 illustrates example performance of an attack generator with a ramp attack policy after each epoch.

FIG. 18 illustrates example performance of an attack generator with a sin attack policy after each epoch.

FIG. 19 illustrates example performance of an attack generator with a pulse attack policy after each epoch.

FIG. 20 illustrate example performance of generated attacks from trained attack generators with different attack policies, according to a study of an example implementation.

FIG. 21 illustrates a schematic description of a networked cyber-physical gas pipeline system.

FIG. 22 illustrates example gas pipeline topologies tested in a study of an example implementation of the present disclosure.

FIG. 23 illustrates a study of the performance of an attack generator with a ramp attack policy applied to a linear topological pipeline system in two epochs.

FIG. 24 illustrates a study of the performance of an attack generator with a ramp attack policy applied to a tree topological pipeline system in two epochs.

FIG. 25 illustrates a study of the performance of an attack generator with a ramp attack policy on the cyclic topological pipeline system in two epochs.

FIG. 26 illustrates the performance of a sample generative attack on a linear topological pipeline system.

FIG. 27 illustrates the performance of a sample generative attack on a tree topological pipeline system.

FIG. 28 illustrates the performance of a sample generative attack on a cyclic topological pipeline system.

FIG. 29 illustrates a table of results for a study of an example implementation of the present disclosure.

FIG. 30 shows the testing scores of the trained generators after batches in the first 2 epochs, according to a study of an example implementation of the present disclosure.

DETAILED DESCRIPTION

Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art. Methods and materials similar or equivalent to those described herein can be used in the practice or testing of the present disclosure. As used in the specification, and in the appended claims, the singular forms “a,” “an,” “the” include plural referents unless the context clearly dictates otherwise. The term “comprising” and variations thereof as used herein is used synonymously with the term “including” and variations thereof and are open, non-limiting terms. The terms “optional” or “optionally” used herein mean that the subsequently described feature, event or circumstance may or may not occur, and that the description includes instances where said feature, event or circumstance occurs and instances where it does not. Ranges may be expressed herein as from “about” one particular value, and/or to “about” another particular value. When such a range is expressed, an aspect includes from the one particular value and/or to the other particular value. Similarly, when values are expressed as approximations, by use of the antecedent “about,” it will be understood that the particular value forms another aspect. It will be further understood that the endpoints of each of the ranges are significant both in relation to the other endpoint, and independently of the other endpoint. While implementations will be described for improving and/or analyzing vulnerabilities in industrial SCADA systems, it will become evident to those skilled in the art that the implementations are not limited thereto, but are applicable for improving the security of communications and control systems in general.

As used herein, the terms “about” or “approximately” when referring to a measurable value such as an amount, a percentage, and the like, is meant to encompass variations of ±20%, ±10%, ±5%, or ±1% from the measurable value.

The term “artificial intelligence” is defined herein to include any technique that enables one or more computing devices or comping systems (i.e., a machine) to mimic human intelligence. Artificial intelligence (AI) includes, but is not limited to, knowledge bases, machine learning, representation learning, and deep learning. The term “machine learning” is defined herein to be a subset of AI that enables a machine to acquire knowledge by extracting patterns from raw data. Machine learning techniques include, but are not limited to, logistic regression, support vector machines (SVMs), decision trees, Naïve Bayes classifiers, and artificial neural networks. The term “representation learning” is defined herein to be a subset of machine learning that enables a machine to automatically discover representations needed for feature detection, prediction, or classification from raw data. Representation learning techniques include, but are not limited to, autoencoders. The term “deep learning” is defined herein to be a subset of machine learning that that enables a machine to automatically discover representations needed for feature detection, prediction, classification, etc. using layers of processing. Deep learning techniques include, but are not limited to, artificial neural network or multilayer perceptron (MLP).

Machine learning models include supervised, semi-supervised, and unsupervised learning models. In a supervised learning model, the model learns a function that maps an input (also known as feature or features) to an output (also known as target or targets) during training with a labeled data set (or dataset). In an unsupervised learning model, the model learns patterns (e.g., structure, distribution, etc.) within an unlabeled data set. In a semi-supervised model, the model learns a function that maps an input (also known as feature or features) to an output (also known as target or target) during training with both labeled and unlabeled data.

Deep learning models, including LLMs, may include artificial neural networks. An artificial neural network (ANN) is a computing system including a plurality of interconnected neurons (e.g., also referred to as “nodes”). This disclosure contemplates that the nodes can be implemented using a computing device (e.g., a processing unit and memory as described herein). The nodes can be arranged in a plurality of layers such as input layer, output layer, and optionally one or more hidden layers. An ANN having hidden layers can be referred to as deep neural network or multilayer perceptron (MLP). Each node is connected to one or more other nodes in the ANN. For example, each layer is made of a plurality of nodes, where each node is connected to all nodes in the previous layer. The nodes in a given layer are not interconnected with one another, i.e., the nodes in a given layer function independently of one another. As used herein, nodes in the input layer receive data from outside of the ANN, nodes in the hidden layer(s) modify the data between the input and output layers, and nodes in the output layer provide the results. Each node is configured to receive an input, implement an activation function (e.g., binary step, linear, sigmoid, tanH, or rectified linear unit (ReLU) function), and provide an output in accordance with the activation function. Additionally, each node is associated with a respective weight. ANNs are trained with a dataset to maximize or minimize an objective function. In some implementations, the objective function is a cost function, which is a measure of the ANN's performance (e.g., error such as L1 or L2 loss) during training, and the training algorithm tunes the node weights and/or bias to minimize the cost function. This disclosure contemplates that any algorithm that finds the maximum or minimum of the objective function can be used for training the ANN. Training algorithms for ANNs include, but are not limited to, backpropagation.

Implementations of the present disclosure include improvements to the security of networked pipeline systems (NPS) and other infrastructure systems. Such systems broadly include “cyber physical systems” that can control the physical elements of a facility or infrastructure. The physical system can be referred to as a “physical plant” and can broadly include machinery, communications systems, networks, and controllers. One example type of infrastructure are networked pipeline systems (NPS) that can play an important role in transporting essential resources such as oil and gas. Increasing integration of information and control techniques in these systems can improve their performance, at the cost of increasing their vulnerability to cyber attacks. Such attacks can have severe consequences on the safety, security, and reliability of pipeline operations. Supervisory control and data acquisition (SCADA) can be used to operate NPS, and provide data collection, state estimation, and decision processes. But SCADA also introduces vulnerabilities. As one example, the state estimation process could be misled maliciously by compromising only a small portion of the IoT-based measurement system [3]. Modifying the control inputs at the automatic control layer can result in catastrophic consequences on physical assets [4].

In order to effectively design and operate security systems, it can be beneficial to model potential attack scenarios and model potential system responses. Such modeling can include generating modeled “attacks,” referred to herein as “attack generation.” Assessing the vulnerability of NPS is a challenging task, as it requires a thorough understanding of the complex interplay between the physical and cyber components of the system and the potential attack scenarios. Existing systems and methods for attack generation and modeling can fail to account for the nonlinear nature of many networked systems (e.g., pipeline systems) and/or require complex or difficult to acquire training data.

Implementations of the present disclosure include improvements to overcome limitations of existing systems for modeling attacks. An example implementation includes an attack generation framework that can integrate physical runtime data of a pipeline system and use a data-driven attack generation approach. The vulnerability analysis can be formulated as determining the presence of feasible attack sets, defined by boundary functions representing the effectiveness and stealthiness of attack signals with respect to the objective and attack detection module. The framework utilizes three data-driven models, including two discriminative models that learn the boundary functions and a generative model that produces elements of the feasible attack set. A loss function defined herein ensures successful attack generation with high probability.

While the implementations described herein are described with reference to industrial and infrastructure systems (e.g., pipelines) it should be understood that implementations of the present disclosure can be used to improve the security and/or operation of any system.

With reference to FIG. 1, an example cybersecurity system is shown according to implementations of the present disclosure. A cybersecurity system 100 can include a trained machine learning model 110 that is configured to assess vulnerabilities of a physical plant 150. The physical plant 150 can include any number of network devices 152 (e.g., computing devices, sensors, controllers, etc.) and any number of physical systems 154 (e.g., pumps, motors, actuators, etc.). The cybersecurity system 100 can use the trained machine learning model 110 to model the physical plant 150 and determine cybersecurity vulnerabilities of the network devices 152 and physical systems 154.

For example, the cybersecurity system 100 can include a controller 112 (e.g., a computing device 1100 as described with reference to FIG. 4). The controller 112 can be configured to receive a trained attack generative model, simulate, by the trained attack generative model, a plurality of attacks on the physical plant, determine an attack effectiveness of each of the plurality of attacks; and control access to the network based on the effectiveness of each of the plurality of attacks.

It should be understood that the physical plant 150 can be any physical system. For example, the physical plant 150 can be a networked pipeline system including pressure sensors, and the cybersecurity system 100 can be configured to secure the pipeline system by controlling access to the network devices 152, and/or securing parts of the physical systems 154 (e.g., securing one or more of the pressure sensors).

As another non-limiting example, the physical plant can include a power grid, and the physical systems 154 can include one or more meters configured to measure parameters of the power grid (e.g., voltage, frequency, current) at nodes in the power grid. Again, the cybersecurity system 100 can be used to model the vulnerabilities of the power grid, and can control access to the network devices 152 (e.g., meters) and/or physical systems 154 of the power grid to improve the security of the power grid.

With reference to FIG. 2A, an example method is shown according to implementations of the present disclosure. The example method can be used to configure the cybersecurity system 100, for example by configuring the trained machine learning model 110 and/or controller 112.

At step 202 the method includes generating a random parameter vector.

At step 204 the method includes generating a random attack dataset.

At step 206 the method include inputting samples of the random attack dataset into a generator. Optionally the generator can be implemented as a deep neural network.

At step 208 the method includes generating a generated attack dataset using the generator. Optionally, the generated attack dataset can include an attack policy that describes the type(s) of attack that can be used. An example category of attack is an attack on one or more sensors of a physical system (e.g., by injecting false sensor data into a network). Non-limiting examples of this type of attack include ramp attacks, where false data is ramped up or down over time.

At step 210 the method includes training a discriminator using the random attack dataset and the generated attack dataset.

At step 212 the method includes training the generator using the trained discriminator and a loss function. An example loss function that can be used is:

$L\left(G\left(z;\theta \right)\right)=\mathrm{ReLU}\left(\underset{z~P_z}{𝔼}\left[D\left(G\left(z;\theta \right)\right)\right]-\left(1-\alpha \right)\right).$

Additional description of loss functions is provided with reference to Example 1, hereto.

FIG. 3A illustrates an example algorithm implementing the method of FIG. 2A. Example 1, herein, further describes the example algorithm of FIG. 3A.

FIG. 2B illustrates another example method according to implementations of the present disclosure. The method of FIG. 2B can be used to configure the cybersecurity system 100, for example by configuring the trained machine learning model 110 and/or controller 112. At step 220 the method includes receiving a generative model parameter vector.

At step 222 the method includes generating sample points from a prior distribution of the generative model parameter vector.

At step 224 the method includes generating attack policy parameters. Non-limiting examples of attack policy parameters include a ramp attack policy, sine attack policy, and pulse attack policy. The present disclosure contemplates that combinations of attack policy parameters can also be used.

At step 226 the method includes simulating simulated attack policy parameters by a system simulation.

At step 228 the method includes training a generative model by the simulated attack policy parameters. Optionally the generative model is a deep neural network, as described in example 2 herein. The deep neural network can optionally select a best loss value by the function: J_best⁽ⁱ⁺¹⁾=min{J_best⁽ⁱ⁾,J(θ_g⁽ⁱ⁺¹⁾)}.

FIG. 3B illustrates an example algorithm implementing the method of FIG. 2B.

It should be appreciated that the logical operations described herein with respect to the various figures may be implemented (1) as a sequence of computer implemented acts or program modules (i.e., software) running on a computing device (e.g., the computing device described in FIG. 4), (2) as interconnected machine logic circuits or circuit modules (i.e., hardware) within the computing device and/or (3) a combination of software and hardware of the computing device. Thus, the logical operations discussed herein are not limited to any specific combination of hardware and software. The implementation is a matter of choice dependent on the performance and other requirements of the computing device. Accordingly, the logical operations described herein are referred to variously as operations, structural devices, acts, or modules. These operations, structural devices, acts and modules may be implemented in software, in firmware, in special purpose digital logic, and any combination thereof. It should also be appreciated that more or fewer operations may be performed than shown in the figures and described herein. These operations may also be performed in a different order than those described herein.

Referring to FIG. 4, an example computing device 1100 upon which the methods described herein may be implemented is illustrated. It should be understood that the example computing device 1100 is only one example of a suitable computing environment upon which the methods described herein may be implemented. Optionally, the computing device 1100 can be a well-known computing system including, but not limited to, personal computers, servers, handheld or laptop devices, multiprocessor systems, microprocessor-based systems, network personal computers (PCs), minicomputers, mainframe computers, embedded systems, and/or distributed computing environments including a plurality of any of the above systems or devices. Distributed computing environments enable remote computing devices, which are connected to a communication network or other data transmission medium, to perform various tasks. In the distributed computing environment, the program modules, applications, and other data may be stored on local and/or remote computer storage media.

In its most basic configuration, computing device 1100 typically includes at least one processing unit 1106 and system memory 1104. Depending on the exact configuration and type of computing device, system memory 1104 may be volatile (such as random access memory (RAM)), non-volatile (such as read-only memory (ROM), flash memory, etc.), or some combination of the two. This most basic configuration is illustrated in FIG. 4 by dashed line 1102. The processing unit 1106 may be a standard programmable processor that performs arithmetic and logic operations necessary for operation of the computing device 1100. The computing device 1100 may also include a bus or other communication mechanism for communicating information among various components of the computing device 1100.

Computing device 1100 may have additional features/functionality. For example, computing device 1100 may include additional storage such as removable storage 1108 and non-removable storage 1110 including, but not limited to, magnetic or optical disks or tapes. Computing device 1100 may also contain network connection(s) 1116 that allow the device to communicate with other devices. Computing device 1100 may also have input device(s) 1114 such as a keyboard, mouse, touch screen, etc. Output device(s) 1112 such as a display, speakers, printer, etc. may also be included. The additional devices may be connected to the bus in order to facilitate communication of data among the components of the computing device 1100. All these devices are well known in the art and need not be discussed at length here.

The processing unit 1106 may be configured to execute program code encoded in tangible, computer-readable media. Tangible, computer-readable media refers to any media that is capable of providing data that causes the computing device 1100 (i.e., a machine) to operate in a particular fashion. Various computer-readable media may be utilized to provide instructions to the processing unit 1106 for execution. Example tangible, computer-readable media may include, but is not limited to, volatile media, non-volatile media, removable media and non-removable media implemented in any method or technology for storage of information such as computer readable instructions, data structures, program modules or other data. System memory 1104, removable storage 1108, and non-removable storage 1110 are all examples of tangible, computer storage media. Example tangible, computer-readable recording media include, but are not limited to, an integrated circuit (e.g., field-programmable gate array or application-specific IC), a hard disk, an optical disk, a magneto-optical disk, a floppy disk, a magnetic tape, a holographic storage medium, a solid-state device, RAM, ROM, electrically erasable program read-only memory (EEPROM), flash memory or other memory technology, CD-ROM, digital versatile disks (DVD) or other optical storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices.

In an example implementation, the processing unit 1106 may execute program code stored in the system memory 1104. For example, the bus may carry data to the system memory 1104, from which the processing unit 1106 receives and executes instructions. The data received by the system memory 1104 may optionally be stored on the removable storage 1108 or the non-removable storage 1110 before or after execution by the processing unit 1106.

It should be understood that the various techniques described herein may be implemented in connection with hardware or software or, where appropriate, with a combination thereof. Thus, the methods and apparatuses of the presently disclosed subject matter, or certain aspects or portions thereof, may take the form of program code (i.e., instructions) embodied in tangible media, such as floppy diskettes, CD-ROMs, hard drives, or any other machine-readable storage medium wherein, when the program code is loaded into and executed by a machine, such as a computing device, the machine becomes an apparatus for practicing the presently disclosed subject matter. In the case of program code execution on programmable computers, the computing device generally includes a processor, a storage medium readable by the processor (including volatile and non-volatile memory and/or storage elements), at least one input device, and at least one output device. One or more programs may implement or utilize the processes described in connection with the presently disclosed subject matter, e.g., through the use of an application programming interface (API), reusable controls, or the like. Such programs may be implemented in a high level procedural or object-oriented programming language to communicate with a computer system. However, the program(s) can be implemented in assembly or machine language, if desired. In any case, the language may be a compiled or interpreted language and it may be combined with hardware implementations.

EXAMPLES

The following examples are put forth so as to provide those of ordinary skill in the art with a complete disclosure and description of how the compounds, compositions, articles, devices and/or methods claimed herein are made and evaluated, and are intended to be purely exemplary and are not intended to limit the disclosure. Efforts have been made to ensure accuracy with respect to numbers (e.g., amounts, temperature, etc.), but some errors and deviations should be accounted for. Unless indicated otherwise, parts are parts by weight, temperature is in ° C. or is at ambient temperature, and pressure is at or near atmospheric.

Example 1

An example implementation of the present disclosure includes an attack generation framework for evaluating the vulnerability of nonlinear networked pipeline systems.

As used in the present example, ⁿ, ₊ denotes the space of real vectors of length n and positive real numbers respectively. The probability triple is denoted by (Ω, , P_Z) [9], where Ω is a sample space containing the set of all possible outcomes, is an event space, and P_z is the associated probability functions of the events in . The study used the symbol to denote the expected value operator. An undirected graph (V, ε) contains vertices, denoted by ε={v₁, v₂, . . . , v_n}, and edges, denoted by ε∈V×V.(v₁, v₂) and (v₁, v₂)∈ε represents an edge in . Consequently, the adjacency matrix, denoted by A(), is a square matrix of size |V|, defined as [10]

$A_{\mathrm{ij}}=\{\begin{array}{cc}1& \mathrm{if}\left(v_i,v_j\right)\in \epsilon \\ 0& \mathrm{otherwise}\end{array}$

An example implementation was studied with a model of a networked pipeline system, including the physical plants, network dynamics, a supervisory control and data acquisition (SCADA) system, and an attack model. FIG. 5 illustrates the system, featuring two control layers. The supervised control layer collects measurement data from IoT sensors across the network, estimates system states, identifies measurement errors, and generates the reference for the local controllers. The automatic control layer includes lower-level controllers of the compressor stations.

A nonlinear creep flow model was used to describe the gas flow in a pipeline segment [11]:

$\begin{array}{cc}\frac{\partial p}{\partial t}=-\frac{c^2\rho }{G}·\frac{\partial Q_n}{\partial x},\frac{\partial p}{\partial x}=-\frac{2f{\rho }^2{c}^2{Q}_n^2}{{\mathrm{DG}}^{2}p}& \left(1\right)\end{array}$

• • where p is the pressure in the pipeline, x is the position of flow, c is the speed of sound in gas [m/s], ρ is the average gas density over the cross-section area of the pipeline [kg/m³], Q_n is the volumetric flow at standard conditions, f is the friction factor, D is the diameter of the pipeline and G=π(D/2)². While the partial differential equation (PDE) model provides an accurate description of the gas flow dynamics for infrequent online optimization, an ordinary differential equation (ODE) model is sufficient to describe the pressure dynamics at nodes [11]. The dynamical pressure model at node i is given as

$\begin{array}{cc}\begin{array}{c}{p}_i=\frac{c^2}{V_i}{\sum }_{j\in {\mathrm{𝒩𝒩}}_i}\sqrt{\frac{❘p_j^2-p_i^2❘}{K_{\mathrm{ij}}}}\mathrm{sgn}\left(p_j-p_i\right)-w_i\\ \mathrm{where}{K}_{\mathrm{ij}}=\frac{64f_{\mathrm{ij}}{c}^2\Delta x_{\mathrm{ij}}}{{\pi }^2{D}_{\mathrm{ij}}^5},V_i=\frac{\pi }{8}{\sum }_{j\in {𝒩}_i}{D}_{\mathrm{ij}}^2\Delta x_{\mathrm{ij}},\end{array}& \left(2\right)\end{array}$

and w_i is the mass flow at pipeline i, Δx_ij denotes the length of pipeline between node i and node j, denote the set of neighborhood pipelines of the pipeline i. Let

$\psi \left(p_i,p_j\right)\stackrel{△}{=}\sqrt{\frac{❘p_j^2-p_i^2❘}{K_{ij}}}\mathrm{sgn}\left(p_j-p_i\right)$ $\mathrm{and}$ $\Psi \left(p\right)\stackrel{△}{=}\left[\begin{array}{cccc}0& \psi \left(p_1,p_2\right)& \dots & \psi \left(p_1,p_n\right)\\ \psi \left(p_2,p_1\right)& 0& \dots & \psi \left(p_2,p_n\right)\\ ⋮& ⋮& \ddots & ⋮\\ \psi \left(p_n,p_1\right)& \psi \left(p_n,p_2\right)& \dots & 0\end{array}\right]$

• • then the model of the entire pipeline network is given by:

$\begin{array}{c}\stackrel{.}{p}=c^2{V}^{-1}\left(A\odot \Psi \left(p\right)\right)1-w\\ =c^2{V}^{-1}\left(A\odot \Psi \left(p\right)\right)1+B_{\sup }{w}_{\sup }-B_{dem}{w}_{dem}\\ \stackrel{△}{=}f\left(p,w\right)\end{array}$

• • where V⁻¹=diag([V₁⁻¹, V₂⁻¹, . . . , V_n⁻¹]), w_sup is a vector of mass inflows at the supply points and w_dem is a vector of mass outflows at the demand points. Clearly w=−B_supw_sup+B_demw_dem, where B_dem and B_sup are appropriately dimensioned matrices such that B_sup^TB_dem=0 and P[−B_sup|B_dem]=I for some projection matrix P. In addition, a dynamical model of the compressor between two nodes i and j is given by [12]:

${\dot{\omega }}_{ij}=\frac{1}{J_{0_i}}\left(u_{ij}+\rho r_2^2{q}_{ij}{\omega }_{ij}\right)$ ${\dot{w}}_{ij}=q_{ij}$

• • where

$\varphi \left(\omega ,q\right)={\left(1+\frac{1}{T_{in}{c}_{p}1000}\left(\rho r_2^2{\omega }^2+k_f{q}^2-\frac{r_1^2}{2}{\left(\omega -\frac{q}{A_c{r}_1{\rho }_{in}}\right)}^2\right)\right)}^{\frac{\gamma }{\gamma -1}},$

ω_ij is the compressor rotor angular velocity [rad/s], q_ij is the mass power flow [kg/s], ϕ(ω, q) is the compressor characterization, u_ij is the mechanical torque [N-m] applied to the rotor shaft or inertia J_0i [kg-m²]. p_i and p_j are the input and output pressures of the compressor, respectively. For the system considered, the compressors are at the extremes of the network. Thus, at each compressor node, one of p_i and p_j will be the ambient pressure. p_i=P_amb for upstream compressors, while p_j=P_amb for downstream compressors.

A supervised control framework was employed to stabilize pipeline pressure and regulate mass flow, compromising a supervised control layer and an automatic control layer. The supervised control layer contains an estimator, a bad data detector (BDD), and a supervised controller. The pipeline system is equipped with hundreds of smart pressure and flow meters, the goal of the estimator at the supervised control layer is to fuse the meter readings and estimate the states of interest p. The measurement model is given as

$\begin{array}{cc}y=g\left(p,w\right)& \left(5\right)\end{array}$

Then an unscented Kalman Filter (UKF) is utilized to perform sensor fusion and state estimation. Firstly, following standard unscented transformation, the study used 2n+1 sigma points to approximate the state p with assumed mean p and covariance P_p as follows:

${\chi }_0=\overline{p},{\chi }_i=\overline{p}+{\left(\sqrt{\left(\lambda +n\right)P_p}\right)}_{i^{\prime }},i=1,\dots ,n$ ${\chi }_{i+n}=\overline{p}+{\left(\sqrt{\left(\lambda +n\right)P_p}\right)}_{i-n},i=n+1,\dots ,2n$

The corresponding weights for the sigma points are given as W₀^m=λ/(n+λ), W₀^c=W₀^m+(1−α²+β), W_i=½(L+λ), and λ=α²(n+κ)−n represents how far the sigma points are away from the state, κ≥0, α∈(0,1], and β=2 is the optimal choice for Gaussian distribution. Assuming p_k−1˜(p_k−1, P_p,k−1), the prediction step is given by:

${\hat{p}}_k^{-}={\sum }_{i=0}^{2n}{W}_{i}g\left(x_{k-1},w_k\right),{\hat{y}}_k^{-}={\sum }_{i=0}^{2n}{W}_i{k,i}_{}$ ${\hat{P}}_{p,k}={\sum }_{i=0}^{2n}{W}_i\left(g\left(X_{k-1},w_k\right)-{\hat{x}}_k\right){\left(g\left(X_{k-1},w_k\right)-{\hat{x}}_k\right)}^T+R$

• • where =f(g(X_k−1, w_k)). The correction step is given by

$\begin{array}{cc}{\hat{p}}_k={\hat{p}}_k^{-}+K_k\left(y_{k,}-{\hat{y}}_{k,}\right),P_{p,k}={\hat{P}}_{p,k}-K_k{\hat{P}}_{y,k}{K}_k^T& \left(6\right)\end{array}$

• • where the Kalman gain is K_k={circumflex over (P)}_py{circumflex over (P)}_y,k⁻¹ with:

${\hat{P}}_{y,k}={\sum }_{i=0}^{2n}{W}_i\left({k,i}_{}-{\hat{y}}_{k,}\right){\left({k,i}_{}-{\hat{y}}_{k,}\right)}^T+Q,$ ${\hat{P}}_{py}={\sum }_{i=0}^{2n}{W}_i\left(X_{k,i}^{\star }-{\hat{p}}_k\right){\left({\left(k,i\right),}_{}-{\hat{y}}_{k,}\right)}^T,$

• • and Q and R are the measurement and process noise covariance matrices respectively. Bad data detector (BDD) is defined as a function D: ⁿ×^p×^m→₊ mapping from the state estimates to a detection residual (i.e. 1-norm or 2-norm residual-based detectors [13], [14]) or a detection likelihood (i.e. χ² detector [15], [8]).

$\begin{array}{cc}{r}_i=D\left({\hat{p}}_i,y_i,w_i\right)& \left(7\right)\end{array}$

Linearizing the models around the equilibrium point (p_eq, w_eq) and discretizing with a fixed time step T_s yields:

$\begin{array}{cc}\Delta p_{k+1}=A\Delta p_k+B\Delta w_k& \left(8\right)\end{array}$

• • where Δp_k=p_k−p_eq, Δw_k=w_k−w_eq, and

${{{A=\frac{\partial f}{\partial p}❘}_{\left(p_{\mathrm{eq}},w_{\mathrm{eq}}\right)}{T}_S+I_n,B=\frac{\partial f}{\partial w}❘}_{\left(p_{\mathrm{eq}},w_{eq}\right)}{T}_S,C=\frac{\partial g}{\partial p}❘}_{P_{\mathrm{eq}}}$

Then, a model predictive controller of horizon length h is utilized to generate the reference mass flow w^r:

${\mathrm{Minimize}}_{\Delta p_{\left[k,k+h\right]},\Delta w_{\left[k,k+h\right]}}:{\sum }_{t=0}^h\frac{1}{2}\Delta p_t^{\top }{M}_1\Delta p_t+{\sum }_{t=0}^{h-1}\frac{1}{2}\Delta w_t^{\top }{M}_2\Delta w_t$ $\mathrm{Subject}\mathrm{to}:\Delta p_{t+1}=A\Delta p_t+B\Delta w_t,t=0,\dots ,h-1$

Thus, the reference mass flow is given as w_k^r=Δw*_[k]+w_eq. Next, at the automatic control layer, a PID controller is used to control the compressors to track the reference mass flow.

$\begin{array}{cc}{u}_{ij}=K_{P}e\left(t\right)+K_I{\int }_0^{t}e\left(\tau \right)d\tau +K_D\frac{\mathrm{de}\left(t\right)}{\mathrm{dt}}& \left(10\right)\end{array}$

• • where e(t)=w_ij−w_ij^r, and the PID gains K_P, K_I, K_D are designed to achieve lim_t→∞∥e(t)∥=0.

In the example study, attacks can be injected through the sensing process and actuation process, and it can be modeled as [15], [14]

$\dot{p}=f\left(p,w+e^u\right)$

To characterize the effect of attacks in the system (11), effectiveness and stealthiness are two common criteria [3]. The study used l₁: ^m×^p→, l₂: ^m×^p→ to denote functions evaluating the effectiveness and stealthiness of the attacks respectively. Then a feasible set of attacks S is defined with given thresholds of effectiveness and stealthiness τ_E, τ_S:

$\begin{array}{cc}=\left\{e^u\in {\Sigma }_{k_u},e^y\in {\sum }_{k_y}❘l_1\left(e^u,e^y\right)\ge {\tau }_E,l_2\left(e^u,e^y\right)\le {\tau }_S\right\}& \left(12\right)\end{array}$

Most research uses the estimation error to represent the effectiveness of sensor attacks such that l₁=∥{circumflex over (x)}_i−x_i∥[13], [14]. The stealthiness function often employs the BDD function l₂=D, where D is defined in (7).

Consequently, the attack generation problem is finding a generative model of the form

$\begin{array}{cc}G\left(z,\theta \right):\left(\Omega ,,P_z\right)\times {ℝ}^{n_g}\to {ℝ}^n& \left(13\right)\end{array}$

• • with a constant tunable parameter vector θ∈^ng, and a prior probability sample space (Ω, , P_Z) with random variables z˜P_z, such that

$\begin{array}{cc}\mathrm{Pr}\left\{G\left(z,\theta \right)\in \right\}\ge \alpha & \left(14\right)\end{array}$

• • for some α∈(0,1). Essentially, G(z, θ) is a generative model for the set if G(z, θ)∈ with a high probability given by the lower bound α. However, obtaining the close-form expressions of l₁ and l₂ can be challenging due to their complex composition of the models from (2) to (10). Moreover, searching for feasible attacks in is at least an NP-hard problem due to the nonlinearity and nonconvexity of the problem. To address these challenges, implementations of the present disclosure can include a data-driven approach that can include two discriminators learning l₁ and l₂ from runtime data, and a generator learning how to solve this NP-hard attack generation problem.

The study introduced a solution to search for feasible attacks in the set defined in (12) using only the runtime data of NPS. The data-driven attack generation framework is illustrated in FIG. 6. The runtime data is obtained from the physical experiment and is used to guide the training of the discriminators. The trained discriminators then supervise the learning process of the generator.

The generator G(z, θ) can be a deep neural network learning the distribution of feasible attacks from sampling noise. It can optionally be trained with the loss function

$\begin{array}{cc}L\left(G\left(z;\theta \right)\right)=\mathrm{ReLU}\left(\underset{z~P_z}{𝔼}\left[D\left(G\left(z;\theta \right)\right)\right]-\left(1-\alpha \right)\right)& \left(15\right)\end{array}$

• • where D function is given by:

$D\left(e\right)=\exp [{\mathrm{Leaky}\mathrm{ReLU}}_s\left({\tau }_E-l_1\left(e\right)\right)$

• • where e=[e^uT e^yT]^T. Then the following minimization program is employed to train the generator and solved by gradient descent.

$\begin{array}{cc}{\theta }^{★}=\underset{\theta }{\arg \min }\underset{z~P_z}{𝔼}L\left(G\left(z;\theta \right)\right)& \left(17\right)\end{array}$

Successful training of the generator can require knowledge of l₁ and l₂ functions. However, it is hard to obtain their closed-form expression. Data-driven approximation offers a practical solution, as it only requires runtime data instead of a high-fidelity model and provides easy calculation of gradients. The study used two deep regression neural networks to approximate l₁ and l₂:

$\begin{array}{cc}{a}_E=f_1\left(e;{\theta }_E\right),a_S=f_2\left(e;{\theta }_S\right)& \left(18\right)\end{array}$

They are both trained with the mean-square-error (MSE) loss function given N runtime effectiveness metric data l₁(e) with the corresponding injection of attacks e:

$\begin{array}{cc}\underset{{\theta }_E}{\mathrm{Minimize}}:\frac{1}{N}{\sum }_{j=1}^N{{l_1\left(e\right)}^{\left(j\right)}-{f_1\left(e;{\theta }_E\right)}^{\left(j\right)}}_2^2& \left(19\right)\end{array}$

• • and given N runtime stealthiness metric data l₂(e) with the corresponding injection of attacks e:

$\begin{array}{cc}\underset{{\theta }_S}{\mathrm{Minimize}}:\frac{1}{N}{\sum }_{j=1}^N{{l_2\left(e\right)}^{\left(j\right)}-{f_2\left(e;{\theta }_S\right)}^{\left(j\right)}}_2^2& \left(20\right)\end{array}$

The generator maps random samples z˜P_z to the parameter vector I∈^p of a pre-defined attack policy π(I). Given I∈^p, the attack policy π(I) is a deterministic time sequence of the injected attacks given by

$\pi \left(I\right)\stackrel{\Delta }{=}\left\{t_0\left(I\right),T\left(I\right),g\left(I\right)\right\}$

It can include the start time of the attack injection t₀(I), the duration of the attack injection T(I), and the attack profile g(I): [t₀(I)t₀(I)+T(I)]→ⁿ. The specific attack policy is assumed to be predetermined for the purpose of the example implementation studied.

Biased learning can be a concern in deep generative models [16], particularly when trained with multiple interrelated data-driven models. The discriminators might learn biased l₁ and l₂ functions based on a limited space of attacks covered by the generative model. This in turn affects the training of the generator network. To improve and address biased learning issues, implementations of the present disclosure can incorporate random and generated attack datasets into the discriminator training process. An example training procedure can include: Receiving Hyperparameters: τ_E, τ_S, α, π, s;

• • Generating a random parameter vector I, implement the attack policy π(I) in system experiment and obtain a random attack dataset {[l₁, l₂, l₃], [a_E, a_S]}_rand;
  • For i in N_epoch:
  • Passing a batch of m samples {z⁽¹⁾, z⁽²⁾, . . . , z^(m)}˜P_z through the generator to obtain [I₁, I₂, I₃]=G(z, θ_i−1);
  • Performing a corresponding attack policy π(G(z, θ_i−1)) in system experiment and obtain a generated attack dataset {[I₁, I₂, I₃], [a_E, a_S]}_gen
  • Training discriminators using both random attack dataset and the current generated attack dataset→(19), 20p;
  • Training the generator using the trained discriminators in the loss function (15), where l₁(e)=f₁(e, θ_Ei) and l₂(e)=f₂(e, θ_Si).

A simulation was performed to evaluate the proposed attack generative model on a 4-node pipeline system, as depicted in FIG. 7. The objective is to regulate the node pressures at their corresponding equilibrium point p_eq, while meeting the demanding mass flow w_dem at node 4 and subject to the flow supply w_h1 and w_h2 at the two wellheads. To this end, the supervised MPC controller (9) generates the reference mass flow w₁ at node 1 using the UKF estimation (6) of the node pressure and the mass flow in the other pipelines. Subsequently, the local PID controller 10) controls the compressor to regulate the pipeline flow supplied to the node 1.

The topology of the network can be represented by the adjacent matrix:

$\begin{array}{cc}A=\left[\begin{array}{cccc}0& 1& 0& 0\\ 1& 0& 1& 0\\ 0& 1& 0& 1\\ 0& 0& 1& 0\end{array}\right]& \left(21\right)\end{array}$

The pipeline model parameters were chosen as c=330 m/s, D₁₂=0.8 m, D₂₃=0.5 m, D₃₄=1.5 m, Δx₁₂=Δx₂₃=Δx₃₄=10 m, f₁₂=f₂₃=f₃₄=0.0025. The compressor model parameters are chosen as J₀=47.7465, ρr₂²=0.0951, T_in c_p=13322.3, k_f=0.05, r₁=20, A_cr₁ρ_in=0.5834, γ=1.2, and A_c/r_c=0.0146

In the simulation, the study targeted the equilibrium point p_eq=[199 194 54 50]^T, corresponding to the equilibrium point of the mass flow w_eq=[−20 −5 −5 13]^T. FIG. 8 illustrates an example of the nominal control performance, where a disturbance d=0.1 is injected into pressure measurements during [100 s, 120 s]. The supervisory control framework successfully regulates the node pressures back to their desired equilibrium point. The UKF is used to estimate 4 node pressure p₁˜p₄ and 3 mass flow w_h1, w_h2, w_dem from their noisy measurements. The BDD can then be defined as

$\begin{array}{cc}r={y-g\left(\hat{x}\right)}_2& \left(22\right)\end{array}$

• • where {circumflex over (x)}=[{circumflex over (p)}^T w_h1 w_h2 w_dem]^T. This BDD residual r is used for the stealthiness metric. And the effectiveness metric is given by the control error

$\begin{array}{cc}e={p-p_{\mathrm{eq}}}_2& \left(23\right)\end{array}$

The objective of the attack generation task is to fool UKF into giving biased estimates {circumflex over (x)}^a leading to small r but big e. The feasible attack set is given particularly as:

$\begin{array}{cc}=\left\{e^y❘e\ge 10,r\le 0\text{.2}\right\}& \left(24\right)\end{array}$

• • where e^y is injected through the measurements as shown in 11). The study used an example attack policy called ramp attack policy π(I)={t₀(I),T(I),g(I)} through each channel, given by:

$\begin{array}{c}{t}_0=t_l+\left(t_r-t_l\right)I_2\\ T=T_l+\left(T_r-T_l\right)I_3\end{array}$

Here ē is the maximum magnitude of the attack signal. The attack injection start time t₀∈[t_l t_r] and the duration time T∈[T_l T_r]. The study use e=0.5, t_l=20 s, t_r=40 s, T_l=0 s, T_r=100 s and the simulation time for generator training is set as 200 s. These exact values of t₀, T, e at runtime are specified by the parameters I∈³ generated by the generator.

Next, the study implemented the example algorithm with τ_E=10, τ_S=0.2, α=0.8, s=0.01 and π given by 25. The generator employed was a deep neural network composed of 4 layers, ReLU, ReLU, Tanh, Sigmoid, with 500, 1000, 500, 3n neurons respectively, where n is the number of nodes under attack. The input size was set as 10. The stealth net includes of 3 ReLU hidden layers and Sigmoid output layer, and the effect net includes of 3 ReLU hidden layers and Linear output layer. The numbers of neurons at layers are 500, 1000, 500, 1 for both effect net and stealth net. They are trained by (20) and 19) respectively. Adams optimizer, in Matlab deep learning toolbox, is used with the learning rate 0.0002, gradient decay factor 0.5, and square gradient decay factor 0.999. The study trained the framework in 5 epochs. In each epoch, the generator trained 10 batches of 5000 batch size, while the discriminators are trained 200 batches of 1000 batch size.

After each epoch of training, the performance of the attack generator is tested in a real system simulation, as shown in FIG. 9. The generated attacks become more stealthy (stealthiness metric value decrease) as the training epoch increases, while their effectiveness is also sacrificed in the first 2 epochs. However, once the stealthiness falls below the threshold 0.2, the generator increases effectiveness from epoch 3 to epoch 5. The F₁ score of generating feasible attacks in the set S defined in 24 is increasing from 0 to 0.86. The training loss curves of discriminators are shown in FIG. 10. It is seen that, at the first 2 epochs, the training of discriminators does not converge yet, so the performance of the attack generator also has bigger deviations as shown in FIG. 9. As the number of epochs increases, the convergence of the generator and two discriminators can arrive simultaneously.

Finally, the study demonstrated the effectiveness of the proposed attack generation method through a time-series performance of generated feasible attacks in FIG. 11. The attacks caused the state estimator to output incorrect estimations of the pipeline flow and current node pressures. Then the controllers were driven to adjust the node pressures towards other malicious equilibrium points. Despite this, the stealthiness value, shown in FIG. 12, remained below the threshold, indicating that the attacks were able to remain undetected.

The example described herein implements a data-driven attack generative model that can be used for NPS vulnerability analysis. The framework includes three interactive data-driven models: two discriminative models to evaluate attack signals and a generative model to generate feasible attacks. The loss function disclosed ensures successful attack generation.

It should be understood that the example in the study is non-limiting, and that other loss functions and/or applications of implementations of the present disclosure are contemplated. For example, the present disclosure contemplates interactive training between the proposed attack generator and a supervised learning-based attack detector, which can lead to a complete automatic exploration of the vulnerability space and the development of a perfect attack detector for the system.

Example 2

A second study was performed on a second example implementation of the present disclosure. The second study included improvements to modeling and design of cyber-physical systems to reduce vulnerability to cyber attacks. As used in the present study, the summation symbol used here is a generalized summation for the appropriate random variable. If the random variable is continuous, then the summation symbol will be considered an integral sign. If the random variable is discrete, then the summation symbol will be the traditional summation of countable quantities. Let x₁, x₂∈. Then, ReLU(x₁)≤ReLU(x₂) implies that x₁≤x₂. Otherwise, x₁≤0 and x₂≤0—in which case ReLU(x₁)=ReLU(x₂)=0.

Modern cyber-physical infrastructure systems (CPIS) are built based on the seamless integration of physical devices and computer systems, facilitating data exchange and control of physical processes [1A]. As the complexity and scale of CPIS increase, network structures are utilized to leverage the resources and capabilities of subsystems, resulting in improved functionality and performance [2A]. In such NCPS, the physical agents, sensors, and actuators are tightly interconnected or coupled through communication networks to gather and analyze data and perform physical actions in response to high-level commands or unwanted external influence. NCPSs are used in various applications, such as smart power grids [3A], industrial pipeline systems [4A], and transportation systems [5A]. The networked and distributed integration of computer systems and physical devices provides numerous benefits, including improved efficiency, reduced costs, and increased safety [6A]. However, networked information transmission presents significant security and privacy challenges as they are vulnerable to cyber-attacks and data breaches [7A]. Therefore, the development of secure and reliable NCPS is essential to ensure their safe and effective operation [8A].

Assessing the vulnerability of NCPS is a challenging task that requires a deep understanding of the complex interplay between the physical and cyber components of the system and the potential attack scenarios. Various techniques to address the vulnerability analysis problem in cyber-physical systems (CPS) exist, but are limited. Such techniques can include mathematical and analytical methods, as well as data-driven and simulation-based approaches. Early researchers incorporated the full system model into the maximization program to generate a feasible attack [9A]-[11A]. In [9A], a false data injection attack (FDIA) was studied against the least-square estimator with a residual-based bad data detector (BDD). The feasibility of FDIA against the Kalman filter with χ{circumflex over ( )}2 detector was studied in [10A]. One study defined vulnerability under sensor attacks as the boundedness of the estimation error and derived sufficient and necessary conditions through analysis of the system's reachable zero dynamics [12A]. A sufficient and necessary condition for insecure estimation under FDIA was also derived for a networked control system in [11A]. To develop more pragmatic attack generation strategies, several constraints are incorporated to capture the attacker's limitations such as limited access to sensors [9A], shallow knowledge of system dynamics [13A], and incomplete knowledge of implemented state estimators [14A]. However, these model-based approaches primarily target linear CPS.

Given the nonlinear nature of NCPSs and the challenges associated with obtaining high-fidelity models, the corresponding vulnerability analysis problem is significantly more complex [15A]. One approach is to analyze vulnerability based on domain knowledge and extensive simulations. For instance, [16A] combined a stochastic adversarial model with a simulation model of interdependent gas-electricity delivery infrastructures to explore operational disruptions caused by cyber-attacks. Another study in [17A] employed hierarchical Limited Stochastic Petri nets and power system network topology to simulate intrusion attack scenarios. Furthermore, a dynamic network security vulnerability assessment method was developed for SCADA systems, taking into account software vulnerabilities and deployed network security defense measures [18A]. However, simulation-based approaches may only be possible for a specific system with specific attack.

Alternatively, data-driven approaches can provide more general solutions. In [19A], [20A], generative adversarial networks (GANs) were trained to learn from existing feasible attack datasets, but their performance relies heavily on the availability of non-generative attack datasets and the representative quality of the training data. Another approach in [15A] utilized data-driven models to approximate the system model, thereby reducing the complexity of solving the optimization-based attack generation problem. However, the optimization-based generator produces a specific attack rather than exploring the vulnerability space.

While existing approaches have been shown to be effective in various applications, they have limitations in capturing the nonlinear nature of NCPSs, relying on exact system model information or prior attack datasets, and being unable to explore vulnerability space. The study disclosed an example implementation including an attack generative framework that overcomes these limitations by integrating knowledge of the underlying physics of NCPS and a data-driven attack generative methodology. A tailored loss function for training the proposed model is given, with backing theoretical guarantees for convergence and probabilistic coverage of the vulnerability space. The example implementation can be applicable to both linear and nonlinear systems without relying on exact model knowledge assumption or prior attack datasets. The example implementation can include a unified attack generation framework for any NCPS. The proposed attack generation system was validated on an IEEE 14-bus system and various gas pipeline systems.

As used in the present example, ⁿ, ₊ and ² denote the Euclidean space of dimension n, set of all positive real numbers, and the Hilbert space of all square summable signals of dimension n respectively. Normal-face lower-case letters (e.g. x∈) are used to represent real scalars, bold-face lower-case letters (e.g. x∈ⁿ) represent vectors, while normal-face upper-case letters (e.g. X∈^m×n) represent matrices. The study used x_i to denote the i th element of the vector x∈ⁿ and x(i) the vector signal x∈ⁿ at time i. dist(θ, Θ) is used to denote the Euclidean distance of a vector θ∈ⁿ to a set Θ*_g⊂ⁿ.

The probability triple is denoted by (Ω, , P_z), where Ω is a sample space containing the set of all possible outcomes, is an event space, and P_z is the associated probability functions of the event z in . The study uses the symbol to denote the expected value operator. The ReLU function is given by

ReLU(x)=max(0,x)

• • and LeakyReLU function is given by

$\mathrm{LeakyReLU}\left(x;s\right)=\max \left(0,x\right)+s\min \left(0,x\right)$

• • for a given parameter s∈(0,1).

(V, ε) denotes an undirected graph with vertices, denoted by ε={v₁, v₂, . . . , v_n}, and edges, denoted by ε∈V×V.(v₁, v₂) and (v₁, v₂)∈ε represents an edge in . Consequently, the adjacency matrix, denoted by A(), is a square matrix of size |V|, defined as [21A]

$A_{\mathrm{ij}}=\{\begin{array}{cc}1& \mathrm{if}\left(v_i,v_j\right)\in ℰ\\ 0& \mathrm{otherwise}\end{array}$

The support of a vector x∈ⁿ is defined as supp(x){i⊆{1, . . . , n}|x_i≠0}. Σ_k{x∈ⁿ∥supp(x)|≤k} denotes the set of k-sparse vectors. The study uses ⊙ to denote element-wise multiplication, i.e.

$\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]\odot \left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]=\left[\begin{array}{cc}{a}_{11}{b}_{11}& a_{12}{b}_{12}\\ a_{21}{b}_{21}& a_{22}{b}_{22}\end{array}\right]$

The study uses ⊗ to denote Kronecker product, i.e. for B∈^m×n,

$\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]\otimes B=\left[\begin{array}{cc}{a}_{11}B& a_{12}B\\ a_{21}B& a_{22}B\end{array}\right]\in {ℝ}^{2m\times 2n}$

The study considers a networked CPS under attack. The networked CPS contains n nodes/agents with communication and computation capability. FIG. 13 illustrates the typical control structure of the networked CPS. The physical plants in the network are equipped with IoT sensors and/or actuators, which are controlled by automatic control units. Communication networks facilitate the sharing of sensing data, while a supervisory control center performs state estimation, bad data detection (BDD), and generates system targets. These targets are then communicated back to the networked plants. However, the transmission of data through the network is susceptible to malicious hackers. Injecting false data can have severe consequences in the physical world.

The example implementation included a system dynamical model. The physical process is modeled by the network dynamics given by:

$\begin{array}{cc}{\stackrel{.}{x}}_i=\sum _{j\in {𝒩}_i}\psi \left(x_i,x_j\right)+B_i{z}_i& \left(1\right)\end{array}$

• • for each node i∈{1, 2, . . . , n}, where x_i∈^m1 is i th node's state vector, z_i∈^m2 is i-th node's actuation state, and denotes the neighborhood of the i-th node/agent, which is the set of nodes physically linked with the i-th node. The physical layer dynamics and the interaction function ψ: ^m1×^m1→^m1 are assumed to satisfy the following properties [22A]:

$1.\psi \left(x_i,x_j\right)=0\mathrm{for}i=j,$ $2.\psi \left(x_i,x_j\right)=-\psi \left(x_j,x_i\right)\mathrm{for}i\ne j,$ $3.\mathrm{The}\mathrm{underlying}\mathrm{graph}\mathrm{is}\mathrm{undirected}.$

Let x=[x₁^T x₂^T . . . x_n^T]^T∈^m1n be an augmented network state vector composing of all nodal state vectors, the study defines the nodal interaction matrix

$\begin{array}{cc}\Psi \left(x\right)\stackrel{△}{=}\left[\begin{array}{cccc}0& \psi \left(x_1,x_2\right)& \dots & \psi \left(x_1,x_n\right)\\ \psi \left(x_2,x_1\right)& 0& \dots & \psi \left(x_2,x_n\right)\\ ⋮& ⋮& \ddots & ⋮\\ \psi \left(x_n,x_1\right)& \psi \left(x_n,x_2\right)& \dots & 0\end{array}\right]& \left(2\right)\end{array}$

Since the underlying graph is undirected, it is shown that the nodal interaction matrix Ψ satisfies the skewsymmetric properties; Ψ(x)^T=−Ψ(x) and η^TΨ(x)η=0 for all X∈^m1n and η∈ⁿ. Consequently, the entire network dynamics are given as

$\begin{array}{cc}\stackrel{.}{x}=\left(A\odot \Psi \left(x\right)\right)1+\mathrm{Bz}& \left(3\right)\end{array}$

• • where A∈^n×n is the adjacency matrix of the underlying undirected graph, B=blkdiag([B₁, B₂, . . . , B_n]) is the block diagonal matrix of all B_i. The network dynamics in (3) is assumed to be controllable, namely; the state-dependent controllability matrix dimension subspace [23A], where L(x) is a state-dependent Laplacian matrix given by

$\begin{array}{cc}{L}_{\mathrm{ij}}\left(x\right)=\{\begin{array}{cc}-\sum _{j\in {𝒩}_i}\frac{\partial \psi \left(x_i,x_j\right)}{\partial x_i}& \mathrm{if}i=j\\ -\frac{\partial \psi \left(x_i,x_j\right)}{\partial x_j}& \mathrm{if}i\ne j\end{array}& \left(4\right)\end{array}$

Let _a denote the index set of actuation nodes, then the actuation state vector z_i=0 if i∉_a, otherwise it follows an actuation dynamics described by the control-affine dynamical model, subject to actuation attacks:

$\begin{array}{cc}{\stackrel{.}{z}}_i=f_i\left(z_i\right)+g_i\left(z_i\right)\left(u_i+e_{u_i}\right)& \left(5\right)\end{array}$

• • where u_i, e_ui∈^p are the i-th node's controlled actuation signal and injected attacks on actuators respectively. f, g are Lipchitz continuous on compact sets.

Next, a measurement model of equipped IoT sensors under attack is given by

$\begin{array}{cc}y=h\left(x\right)+e_y+v& \left(6\right)\end{array}$

• • where y, e_y, v∈^m are the measurement, sensor attacks, and sensor noise respectively. Let =supp(e_u)⊂{1, 2, . . . , pn}, =supp(e_y)⊂{1, 2, . . . , m} denote the attack supports at actuation and sensing processes respectively. The study assumes that the attack signals are sparse. This is typically done in literature [24A] to capture limited resources available to the attacker. Thus, e_u∈Σ_ku and e_y∈Σ_ky, where k_u<p, k_y<m.

The example implementation further included a nominal control design. Here, the study provides a control design for the nominal system without attacks (i.e. e_u=0, e_y=0). The materials in this subsection can be inferred from standard literature on nonlinear control design and are presented for the sake of completeness. The control design comprises 2 layers; (1) the low-level actuator control (Level 1), and (2) the supervisory reference generator (Level 2).

The low-level controllers are closed-loop schemes u_i: ^m2→^p designed to ensure that the actuator state z_i tracks a given reference z_i^r∈^p such that

$\begin{array}{cc}\tilde{z}\left(t\right)\le \tilde{z}\left(0\right)e^{-\lambda t}& \left(7\right)\end{array}$

• • where {tilde over (z)}_i=z_i−z_i^r is the associated tracking error, λ>0 is the convergence rate. As a result, it is assumed that there exists a continuous differentiable, positive definite function V: ^m2→ which satisfies [25A]

$\underset{\_}{\alpha }\left(\tilde{z}\right)\le V\le \stackrel{\_}{\alpha }\left(\tilde{z}\right),\frac{\partial V}{\partial \tilde{z}}\stackrel{.}{\stackrel{~}{z}}\le -\alpha \left(\tilde{z}\right)$

Thus, there exists a stabilizing controller of the form u=K(V({tilde over (z)}))—for example using the Sontag formula [26A]

$\begin{array}{cc}K\left(V\left(\tilde{z}\right)\right)=\{\begin{array}{cc}-\frac{L_{f}V\left(\tilde{z}\right)+\sqrt{{\left[L_{f}V\left(\tilde{z}\right)\right]}^2+{\left[L_{g}V\left(\tilde{z}\right)\right]}^4}}{L_{g}V\left(\tilde{z}\right)}& \mathrm{if}{L}_{g}V\left(\tilde{z}\right)\ne 0\\ 0& \mathrm{if}{L}_{g}V\left(\tilde{z}\right)=0\end{array}\mathrm{where}{L}_{f}V\left(\tilde{z}\right)\stackrel{△}{=}\langle \frac{\partial V}{\partial \tilde{z}},f\left(\tilde{z}\right)\rangle ,L_{g}V\left(\tilde{z}\right)\stackrel{△}{=}\langle \frac{\partial V}{\partial \tilde{z}},g\left(\tilde{z}\right)\rangle & \left(8\right)\end{array}$

are the Lie derivatives of V along f and g respectively.

The supervisory control generates the state reference z^r for the low-level controllers in order to regulate the network states x at the equilibrium point. Consistent with the time-scale separation inherent in a typical NCPS, the study assumes that the lowlevel control layer is much faster than the supervisory layer. Thus, for the purpose of designing a supervisory target generator, the study linearized the system dynamics about the operating point [27A]. The operating point is a known special equilibrium point for which the nominal operation of the CPS is designed. Examples of operating points include the nominal frequency and bus voltage of an AC transmission system (frequency of 60 Hz and voltage ranging from 230 KV to 500 KV for a typical transmission system in the United States). Although the nonlinear dynamics in (3) would have multiple equilibrium points, the study uses the known operating point as the equilibrium point of interest. Let x_eq be the equilibrium point of interest corresponding to the equilibrium points of inputs u_eq. Next, a linear model is used to approximate the network dynamical model in (3) around x_eq:

$\begin{array}{cc}\stackrel{.}{\stackrel{~}{x}}=-L\left(x_{\mathrm{eq}}\right)\tilde{x}+\mathrm{Bz}& \left(9\right)\end{array}$

• • where {tilde over (x)}=x−x_eq is the regulation error and L(x) is defined in (4).

Next, given a sample time step T_s, a discrete form of the linear model in (9) is given as:

$\begin{array}{cc}\tilde{x}\left(k+1\right)=\underset{A_d}{\underbrace{\left(-L\left(x_{\mathrm{eq}}\right)T_s+I\right)}}\tilde{x}\left(k\right)+\underset{B_d}{\underbrace{{\mathrm{BT}}_s}}z\left(k\right).& \left(10\right)\end{array}$

Consequently, a model predictive control (MPC) scheme is utilized to regulate the network states at the equilibrium point x_eq:

$\underset{z_{\left[t,t+h-1\right],}^r}{\mathrm{Minimize}}:{\sum }_{k=0}^h\frac{1}{2}{\tilde{x}\left(k\right)}^T{M}_1\tilde{x}\left(k\right)+{\sum }_{k=0}^{h-1}\frac{1}{2}{z^r\left(k\right)}^T{M}_2{z}^r\left(k\right)$ ${\tilde{x}}_{\left[t,t+h\right]}$ $\mathrm{Subject}\mathrm{to}:\tilde{x}\left(k+1\right)=A_d\tilde{x}\left(k\right)+B_d{z}^r\left(k\right),$ $k=0,\dots ,h-1,$

The goal of the MPC program is to minimize the deviation from the operating point and the corresponding control energy using weight matrices M₁, M₂ respectively. The discrete error dynamics is imposed as an equality constraint. The initial constraint {tilde over (x)}(0)={tilde over (x)}_t and terminal constraint {tilde over (x)}(k)=0 are imposed as another equality constraints. Since the network state x is not measurable, the regulation error is given by {tilde over (x)}(t)={circumflex over (x)}(t)−x_eq, where the estimate of the network states {circumflex over (x)}_t is obtained using the following state estimator, given sensor measurements y.

The state estimator is a minimizer of the difference between model measurements and sensing readings in T horizon, according to the program

$\underset{{\hat{x}}_{\left[t-T+1,t\right]}}{\mathrm{Minimize}}:{\sum }_{k=t-T+1}^t{y\left(k\right)-h\left(\hat{x}\left(k\right)\right)}_2^2$ $\mathrm{Subject}\mathrm{to}:\tilde{x}\left(k+1\right)=A_d\tilde{x}\left(k\right)+B_{d}z\left(k\right)=0$

As the estimator estimates the network states from the measurements, a BDD monitors the state estimation process to detect any false inputs. It is defined as a function D_∈: ^m1n×^m→₊ mapping from the state estimates to a detection residual (i.e. 1-norm or 2-norm residual-based detectors [9A], [28A]) or detection likelihood (i.e. χ² detector [10A], [15A]).

$\begin{array}{cc}{r}_t=D_{ϵ}\left({\hat{x}}_t,y_t,u_t\right)& \left(13\right)\end{array}$

The example implementation further included analyses of a vulnerability analysis problem.

Since obtaining a high-fidelity model for complex networked systems is challenging even for experienced practicing engineers and operators, the system models described in the present example are assumed to be unknown. However, for the purpose of vulnerability analysis, the study assumed that all measurement and actuation channels can be accessed. In other words, the study assumes all actuators and sensors are potential targets for attacks in order to cover all possible vulnerable scenarios. Additionally, the study assumes that the effectiveness and stealthiness of the injected attacks could be measured. Stealthiness refers to the potential to bypass BDD, while effectiveness represents the closeness to the intended degradation of system performance [11A]. Various metrics have can be used to measure effectiveness and stealthiness. For example, estimation error was used as the effectiveness metric while the estimation residual was used as the stealthiness metric in [10A], [14A]. [28A] used a critical measurement, such as the sum of the water levels in a water tank system, as the effectiveness metric. The present example can use any of these metrics, and/or allow the system operator to select metrics and consider them as hyperparameters to the vulnerability analysis problem. Next, the study defines a vulnerable NCPS as follows.

Definition 1 (Vulnerable System). For the NCPS in (3), (5) and (6), let l₁: ²→, l₂: ²→ be the corresponding Lipschitz effectiveness and stealthiness functionals respectively. Consider a vulnerability set

$\left({\tau }_E,{\tau }_S\right)=\left\{e_u\in {\sum }_{k_u},e_y\in {\sum }_{k_y}❘l_1\left(e_u,e_y\right)\ge {\tau }_E,l_2\left(e_u,e_y\right)\le {\tau }_S\right\},$

where τ_E, τ_S are the corresponding effectiveness and stealthiness thresholds. Then

• • 1. The NCPS is said to be vulnerable if (τ_E, τ_S)∉Ø;
  • 2. If (τ_E, τ_S)∉Ø, the attacks e_u, e_y are said to be feasible for given thresholds τ_E and τ_S if {e_u, e_y}∈(τ_E, τ_S)

Specifying the l₁ and l₂ as Lipschitz functionals is justifiable as they are often given by energy-like functions expressed in terms of certain norms and/or inner products. This can be achieved by either employing a linear model approximation [29A] or utilizing data-driven universal approximation [15A].

To assess the vulnerability of the NCPS in (3), (5) and (6), which is equivalent to certifying the nonemptyness of , it is sufficient to generate any feasible e_u, e_y. Therefore, the vulnerability analysis can be formulated as an attack generation problem [30A]. The study shows that the domain of (τ_E, τ_S) is a function space, thereby rendering the attack generation problem infinite dimensional. To circumvent this, the study employs an attack policy, which is a parameterized basislike signal vector used as a heuristic to simplify the search space for the generative model. This essentially transforms the original infinite-dimensional search space into a finitedimensional parameter search. Next, the study formalizes this definition of attack policy.

Definition 2 (Attack Policy). An attack policy π(ϕ) is a mapping from policy parameter vector ϕ∈^nϕ to a deterministic time sequence of the injected attack vectors e.

Thus, this results in a more conservative vulnerability set

$\left({\tau }_E,{\tau }_S\right)=\left\{\varphi \in {ℝ}^{n_{\varphi }}❘l_1\left(\pi \left(\varphi \right)\right)\ge {\tau }_E,l_2\left(\pi \left(\varphi \right)\right)\le {\tau }_S\right\}\underset{\_}{\subset }\left({\tau }_E,{\tau }_S\right)$

• • whose domain is a finite-dimensional Euclidean space. Consequently, to explore the set (τ_E, τ_S), the study developed a generative model for the parameters of the attack policy such that the resulting attack is feasible with a high probability.

Definition 3 (Generative model). A generative model is a mapping of the form

$\begin{array}{cc}G\left(z,\theta \right):\left(\Omega ,,P_z\right)\times {ℝ}^{n_{\theta }}\to {ℝ}^{n_{\varphi }}& \left(15\right)\end{array}$

• • with a tunable parameter vector θ∈^nθ and a probability prior (Ω, , P_z) with random variables z˜P_z, which satisfies

$\begin{array}{cc}\mathrm{Pr}\left\{G\left(z,\theta \right)\in \right\}\ge \eta & \left(16\right)\end{array}$

• • for a given η∈(0,1).Essentially, G(z, θ) is a generative model for the set if G(z, θ)∈ with a high probability given by the lower bound η. The attack policy is ultimately deployed together with a trained generative model. The effectiveness and stealthiness, measured by the functionals l₁ and l₂, are used to train the generative model. Therefore, the goal of training the generative model is to find any feasible parameter vector θ*_g which satisfies (16). The next results gives a sufficient condition for this requirement.

Theorem 1. Given a set whose boundary is of measure zero. If

$\begin{array}{cc}𝔼\left[D\left(G\left(z,\theta \right)\right)\right]\le 1-\eta & \left(17\right)\end{array}$

• • where D: ^nϕ→ is a function which satisfies

$\begin{array}{cc}D\left(\varphi \right)\{\begin{array}{cc}\in \left[0,1\right]& \mathrm{if}\varphi \in \\ >1& \mathrm{otherwise}\end{array}& \left(18\right)\end{array}$

• • then

$\mathrm{Pr}\left\{G\left(z,\theta \right)\in \right\}\ge \eta$

Proof. By definition of the expectation operator, for any random variable x,¹

$\begin{array}{c}𝔼\left[D\left(x\right)\right]=\sum _{X\in {ℝ}^n}D\left(X\right)\mathrm{Pr}\left\{x=X\right\}\\ =\sum _{X\in 𝒮}D\left(X\right)\mathrm{Pr}\left\{x=X\right\}+\sum _{X\notin 𝒮}D\left(X\right)\mathrm{Pr}\left\{x=X\right\}\end{array}$

Since D(x)≥0, it follows that

$𝔼\left[D\left(x\right)\right]\ge \sum _{X\notin 𝒮}D\left(X\right)\mathrm{Pr}\left\{x=X\right\}$

In addition, since D(x)>1 for all x∉, [D(x)] can be further lower bounded as

$𝔼\left[D\left(x\right)\right]\ge \sum _{X\notin 𝒮}\mathrm{Pr}\left\{x=X\right\}=\mathrm{Pr}\left\{x\notin \right\}$

Replacing x with G(z, θ) yields

[D(G(z,θ))]≥Pr{G(z,θ)∉}

Using the inequality in (17), it follows that

$\begin{array}{c}\mathrm{Pr}\left\{G\left(z,\theta \right)\in 𝒮\right\}=1-\mathrm{Pr}\left\{G\left(z,\theta \right)\notin \right\}\\ \ge 1-𝔼\left[D\left(G\left(z,\theta \right)\right)\right]\\ \ge 1-\left(1-\eta \right)\end{array}$ $\ge \eta$

Consequently, for a given threshold η∈(0,1), the study defined the feasible parameter set Θ_g(η):

${\Theta }_g^{*}\left(\eta \right)\stackrel{△}{=}\left\{{\theta }_g\in {ℝ}^{n_g}❘\underset{z~P_z}{𝔼}\left[D\left(G\left(z,{\theta }_g\right)\right)\right]\le 1-\eta \right\}$

By utilizing any function D which satisfies the discrimination property in (18), the nonconvex vulnerability analysis problem can be transformed to a training problem for the generative model G(z, θ). For this, the study considers the example loss function

$\begin{array}{cc}J\left(\theta \right)=\mathrm{ReLU}\left(\underset{z~P_z}{𝔼}\left[D\left(G\left(z;\theta \right)\right)\right]-\left(1-\eta \right)\right)& \left(20\right)\end{array}$

Theorem 2. Given η∈(0,1), if Θ*_g(η) is nonempty, then

$\underset{\theta \in {ℝ}^{n_g}}{\arg \min }J\left(\theta \right)\subseteq {\Theta }_g^{*}\left(\eta \right)$

Proof. Let

${\theta }^{*}\in \underset{\theta \in {ℝ}^{n}g}{\arg \min }J\left(\theta \right).$

Then, given any arbitrary θ₁∈^ng, it follows from the optimality of θ* that

$\mathrm{ReLU}\left(\underset{z~P_z}{𝔼}\left[D\left(G\left(z;{\theta }^{*}\right)\right)\right]-\left(1-\eta \right)\right)\le \mathrm{ReLU}\left(\underset{z~P_z}{𝔼}\left[D\left(G\left(z;{\theta }_1\right)\right)\right]-\left(1-\eta \right)\right)$

Now, due to the order preservation of the ReLU function, either

$\begin{array}{cc}{\theta }_1\in \underset{\theta \in {ℝ}^{n_g}}{\arg \min }J\left(\theta \right)\mathrm{or}\underset{z~P_z}{𝔼}\left[D\left(G\left(z;{\theta }^{*}\right)\right)\right]\le \underset{z~P_z}{𝔼}\left[D\left(G\left(z;{\theta }_1\right)\right)\right]& \left(21\right)\end{array}$

(21) holds for all

${\theta }_1\notin \underset{\theta \in {ℝ}^{n_g}}{\arg \min }J\left(\theta \right).$

Thus, for any

${\theta }_g\in {\Theta }_g^{*}\left(\eta \right)\backslash \underset{\theta \in {ℝ}^{n_g}}{\arg \min }J\left(\theta \right),$

it follows that

$\underset{z~P_z}{𝔼}\left[D\left(G\left(z;{\theta }^{*}\right)\right)\right]\le \underset{z~P_z}{𝔼}\left[D\left(G\left(z;{\theta }_g\right)\right)\right]\le 1-\eta$

• • which implies that θ*∈Θ*_g(θ). Thus far, the study has formulated the vulnerability analysis problem as an attack generation problem solved by training a generative network with a custom loss function. The next development focuses on the training algorithm and associated convergence properties.

Results

The training convergence of the generative model with the loss function in (20), using the specific discrimination function of the form:

D(ϕ)=D(τ_E,τ_S,l₁∘π(ϕ),l₂∘π(ϕ))

• • where the function D is Lipschitz in each argument with Lipschitz constants not bigger than L_D>0. The training process is shown schematically in FIG. 14.

First, the study analyzed the convergence guarantee when l₁, l₂ are known exactly. However, calculating the gradient of such a nonlinear networked dynamical system can be computationally expensive. As a result, the study used universal approximators to learn the functions l₁ and l₂ from runtime data. The study also shows that the convergence property is preserved using the approximated versions {circumflex over (l)}₁ and {circumflex over (l)}₂. To obtain a parameter vector within the feasible set Θ*_g(η), the study employed the iterative scheme:

$\begin{array}{cc}{\theta }_g^{\left(i+1\right)}={\theta }_g^{\left(i\right)}-{\lambda }_i{h}^{\left(i\right)}& \left(23\right)\end{array}$

• • where h⁽ⁱ⁾∈∂J(θ_g⁽ⁱ⁾) is a subgradient of J at θ_g⁽ⁱ⁾ which satisfies the inequality

$\begin{array}{cc}J\left(\theta \right)\ge J\left({\theta }_g^{\left(i\right)}\right)+h^{\left(i\right)\top }\left(\theta -{\theta }_g^{\left(i\right)}\right)& \left(24\right)\end{array}$

• • for all θ. Then, at each iterative step, the study set

$\begin{array}{cc}{J}_{\mathrm{best}}^{\left(i+1\right)}=\min \left\{J_{\mathrm{best}}^{\left(i\right)},J\left({\theta }_g^{\left(i+1\right)}\right)\right\}& \left(25\right)\end{array}$

• • and I_best⁽ⁱ⁺¹⁾=i if J(θ_g⁽ⁱ⁺¹⁾)=J_best⁽ⁱ⁺¹⁾, otherwise I_best⁽ⁱ⁺¹⁾=I_best⁽ⁱ⁾. Accordingly, the study set the parameter estimate as

${\hat{\theta }}_g^{\left(i\right)}={\theta }_g^{\left(I_{\mathrm{best}}^{\left(i\right)}\right)}.$

Proposition 1. Suppose there exists a compact subsets Θ_g⊂^nθ, Φ⊂^nϕ, and real constants L₁>0, L₂>0, L_G>0 such that the generative model G(z, θ) and the functions l₁∘π: ^nϕ→ and l₂∘π: ^nϕ→ satisfy the Lipschitz conditions:

$\begin{array}{cc}\underset{z~P_z}{𝔼}\left[G\left(z,{\theta }_1\right)-G\left(z,{\theta }_2\right)\right]\le L_G{\theta }_1-{\theta }_2& \left(26\right)\end{array}$ $\begin{array}{cc}❘l_1\circ \pi \left({\varphi }_1\right)-l_1\circ \pi \left({\varphi }_2\right)❘\le L_1{\varphi }_1-{\varphi }_2& \left(27\right)\end{array}$ $\begin{array}{cc}❘l_2\circ \pi \left({\varphi }_1\right)-l_2\circ \pi \left({\varphi }_2\right)❘\le L_2{\varphi }_1-{\varphi }_2& \left(28\right)\end{array}$

• • for all θ₁, θ₂∈Θ_g and all ϕ₁, ϕ₂∈Φ. Then, the loss function in (20), with discrimination function D in (22), satisfies the Lipschitz condition

$\begin{array}{cc}❘J\left({\theta }_1\right)-J\left({\theta }_2\right)❘\le L_D\left(L_1+L_2\right)L_G{\theta }_1-{\theta }_2& \left(29\right)\end{array}$

• • for all θ₁, θ₂∈Θ_g.Proof. According to (26),

$\underset{z~P_z}{𝔼}\left[G\left(z,\theta \right)\right]$

is Lipschitz continuous. Thus, for any compact set Θ_g⊂^nθ, the study set Φ as its image under

$\underset{z~P_z}{𝔼}\left[G\left(z,\theta \right)\right].$

Thus, for all θ₁, θ₂∈Θ_g, the following follows from the Lipschitz continuity of the ReLU function

$❘J\left({\theta }_1\right)-J\left({\theta }_2\right)❘\le ❘\underset{z~P_z}{𝔼}\left[D\left(G\left(z,{\theta }_1\right)\right)\right]-\underset{z~P_z}{𝔼}\left[D\left(G\left(z,{\theta }_2\right)\right)\right]❘\le ❘\underset{z~P_z}{𝔼}\left[D\left(G\left(z,{\theta }_1\right)\right)-D\left(G\left(z,{\theta }_2\right)\right)\right]❘$

• • where the last inequality follows from the linearity of the expectation operator. Using the Jensen's inequality, and the conditions in (27) and (28), yields

$❘J\left({\theta }_1\right)-J\left({\theta }_2\right)❘\le \underset{z~P_z}{𝔼}\left[❘D\left(G\left(z,{\theta }_1\right)\right)-D\left(G\left(z,{\theta }_2\right)\right)❘\right]\le \underset{z~P_z}{𝔼}\left[L_D\left(L_1+L_2\right)G\left(z,{\theta }_1\right)-G\left(z,\theta \right)\right]\le L_D\left(L_1+L_2\right)\underset{z~P_z}{𝔼}\left[G\left(z,{\theta }_1\right)-G\left(z,{\theta }_2\right)\right]\le L_D\left(L_1+L_2\right)L_G{\theta }_1-{\theta }_2$

Remark 2. As an example of (22), it is not hard to check that

$D_s\left(\varphi \right)=\exp (\mathrm{LeakyReLU}\left({\tau }_E-l_1\circ \pi \left(\varphi \right);s\right)$

• • satisfies the discrimination property in (18) for any s∈(0, 1). Also, Proposition 1 holds with L_D=e^b, where

$b=\underset{\varphi \in 𝒮\left({\tau }_E,{\tau }_S\right)}{\sup }\left(❘{\tau }_E-l_1\circ \pi \left(\varphi \right)❘+❘{\tau }_s-l_2\circ \pi \left(\varphi \right)❘\right)$

Next, the study presents results on the asymptotic convergence of the scheme in (23).

Theorem 3. Consider a generative model G(z, θ_g) and discriminator D(ϕ) which satisfy the Lipschitz conditions in Proposition 1. Given any initial guess θ_g⁽⁰⁾∈Θ_g. Suppose that the feasible parameter set satisfies the inclusion Θ*_g(η)≤Θ_g, and that the step size sequence of the iterative scheme in (23) is chosen to satisfy the nonsummable square-summable property

$\begin{array}{cc}\sum _{i=0}^{\infty }{\lambda }_i=\infty ,\sum _{i=0}^{\infty }{\lambda }_i^2\le M& \left(31\right)\end{array}$

• • for some real number M<∞. Then, for any ϵ>0,

$\begin{array}{cc}{J}_{\mathrm{best}}^{\left(i\right)}-J\left({\theta }_g^{*}\right)\le ϵ\mathrm{for}\mathrm{all}i\ge N_{ϵ}& \left(32\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}{N}_{ϵ}=\min \left\{i\in \mathbb{Z}❘\sum _{j=0}^i{\lambda }_j\ge \frac{{{\theta }_g^{\left(0\right)}-{\theta }_g^{*}}^2+{L_D^2\left(L_1+L_2\right)}^2{L}_G^{2}M}{2ϵ}\right\}& \left(33\right)\end{array}$ $\mathrm{and}$ ${\theta }_g^{*}\stackrel{△}{=}\underset{\theta \in {\Theta }_g^{*}\left(\eta \right)}{\arg \min }\left\{{\theta }_g^{\left(0\right)}-\theta \right\}$

• • is the projection of the initial guess to the feasible parameter set in (19), L_D and L_G are the Lipschitz constants given in (29).Proof. From (23) and (24),

$\begin{array}{c}{{\theta }_g^{\left(i+1\right)}-{\theta }_g^{*}}^2={{\theta }_g^{\left(i\right)}-{\lambda }_i{h}^{\left(i\right)}-{\theta }_g^{*}}^2\\ ={{\theta }_g^{\left(i\right)}-{\theta }_g^{*}}^2-2{\lambda }_i{h}^{\left(i\right)\top }\left({\theta }_g^{\left(i\right)}-{\theta }_g^{*}\right)+{\lambda }_i^2{h^{\left(i\right)}}^2\le \\ {{\theta }_g^{\left(i\right)}-{\theta }_g^{*}}^2-2{\lambda }_i\left[J\left({\theta }_g^{\left(i\right)}\right)-J\left({\theta }_g^{*}\right)\right]+{\lambda }_i^2{h^{\left(i\right)}}^2\end{array}$

Iterating the last inequality over {0, . . . , i+1} yields

${{\theta }_g^{\left(i+1\right)}-{\theta }_g^{*}}^2\le {{\theta }_g^{\left(0\right)}-{\theta }_g^{*}}^2-2\sum _{j=0}^i{\lambda }_j\left[J\left({\theta }_g^{\left(j\right)}\right)-J\left({\theta }_g^{*}\right)\right]+\sum _{j=0}^i{\lambda }_i^2{h^{\left(j\right)}}^2$

Since ∥θ_g⁽ⁱ⁺¹⁾−θ*_g∥²≥0, it follows that

$\begin{array}{cc}2\sum _{j=0}^i{\lambda }_j\left[J\left({\theta }_g^{\left(j\right)}\right)-J\left({\theta }_g^{*}\right)\right]\le {{\theta }_g^{\left(0\right)}-{\theta }_g^{*}}^2+\sum _{j=0}^i{\lambda }_i^2{h^{\left(j\right)}}^2& \left(34\right)\end{array}$

On the other hand, the following inequalities hold

$\sum _{j=0}^i{\lambda }_j\left[J\left({\theta }_g^{\left(j\right)}\right)-J\left({\theta }_g^{*}\right)\right]\ge \left(\sum _{j=0}^i{\lambda }_j\right)\underset{j=1,\dots ,i}{\min }\left[J\left({\theta }_g^{\left(j\right)}\right)-J\left({\theta }_g^{*}\right)\right]$

• • where (35) follows from (25) by noting that J_best⁽ⁱ⁾=min{J(θ_g⁽⁰⁾), . . . , J(θ_g⁽ⁱ⁾)}. Combining (34) and (35) yields:

$J_{\mathrm{best}}^{\left(i\right)}-J\left({\theta }_g^{*}\right)\le \frac{{{\theta }_g^{\left(0\right)}-{\theta }_g^{*}}^2+{\sum }_{j=0}^i{\lambda }_j^2{h^{\left(j\right)}}^2}{2{\sum }_{j=0}^i{\lambda }_j}$

From Proposition 1, it is seen that the loss function J(θ) is uniformly Lipschitz over Θ*_g(η) with the constant L_D(L₁+L₂)L_G. Thus, the subgradient h⁽ⁱ⁾ can be upper bounded as:

$\begin{array}{cc}{h^{\left(i\right)}}^2\le L_D\left(L_1+L_2\right)L_G& \left(36\right)\end{array}$

leading to the inequality:

$J_{\mathrm{best}}^{\left(i\right)}-J\left({\theta }_g^{*}\right)\le \frac{{{\theta }_g^{\left(0\right)}-{\theta }_g^{*}}^2+{L_D^2\left(L_1+L_2\right)}^2{L}_G^2{\sum }_{j=0}^i{\lambda }_j^2}{2{\sum }_{j=0}^i{\lambda }_j}\le \frac{{{\theta }_g^{\left(0\right)}-{\theta }_g^{*}}^2+{L_D^2\left(L_1+L_2\right)}^2{L}_G^{2}M}{2{\sum }_{j=0}^i{\lambda }_j}\le ϵ\mathrm{for}\mathrm{all}i\ge N_{ϵ}$

Remark 3. Since ϵ is arbitrary and J⁽ⁱ⁾≥J(θ*), it follows that lim_i→∞J_best⁽ⁱ⁾=J(θ*), which according to the feasibility set defined in (19), implies that lim_i→∞θ_g⁽ⁱ⁾∈Θ_g(η) for any given η.

Next, to reduce the computational burden of the gradient evaluation, the study approximated the functions l₁∘π, l₂∘π using universal approximators and develop a corresponding training algorithm for the generative model. The two data-driven models {circumflex over (l)}₁, {circumflex over (l)}₂ are assumed to satisfy the universal approximation properties [31A], [32A]: For any ϵ₁, ϵ₂∈₊, there exists parameter vectors β₁ and β₂ such that

$\underset{\theta \in {\Theta }_g}{\sup }\max \left\{❘{\hat{l}}_1\left(\theta ;{\beta }_1\right)-l_1\circ \pi \left(\theta \right)❘,{\hat{l}}_1-l_1\right\}\le {ϵ}_1$

• • where I_i∈∂(l₁∘π(θ)) and Î_i∈∂{circumflex over (l)}₁(θ; β_i), i∈1,2 are respective subgradients of functions l_i∘π and their estimates {circumflex over (l)}_i(θ; β_i). The inequalities in (37) are expressed (in Sobolev space) in terms of essential supremum norm up to the first weak derivative. This has been used extensively to study the approximation properties of deep neural networks [33A].

Corollary 3.1 (Convergence with Approximated {circumflex over (l)}₁, {circumflex over (l)}₂). Suppose the data-driven models {circumflex over (l)}₁(θ_E⁽ⁱ⁾), {circumflex over (l)}₂(θ_S⁽ⁱ⁾) satisfy the universal approximation properties in (37), where θ_E⁽ⁱ⁾ and θ_S⁽ⁱ⁾ are optimal parameters at epoch i. Then, by following Algorithm 1 with the nonsummable square-summable property in (31), for any ϵ>0,

$J_{best}^{\left(i\right)}-J\left({\theta }_g^{★}\right)\le ϵ\mathrm{for}\mathrm{all}i\ge {\overline{N}}_{ϵ}$ $\mathrm{where}$ ${\overline{N}}_{ϵ}=\min \left\{i\in \mathbb{Z}❘{\sum }_{j=N_{ϵ}+1}^i{\lambda }_j\ge \frac{L_D^2\left(L_1+L_2+\frac{{ϵ}_1+{ϵ}_2}{2d_g}\right)\left(\frac{{ϵ}_1+{ϵ}_2}{d_g}\right)L_G^{2}M}{ϵ}\right\},$

and N_ϵ is the constant given in (33).

Proof. Using the data-driven models {circumflex over (l)}₁ and {circumflex over (l)}₂, the iterative scheme becomes

${\theta }_g^{\left(i+1\right)}={\theta }_g^{\left(i\right)}-{\lambda }_i{\hat{h}}^{\left(i\right)}$

• • where ĥ⁽ⁱ⁾∈∂Ĵ(θ_g⁽ⁱ⁾) is an estimation of the subgradient of the objective function evaluated using the approximate datadriven models. Here, the resulting estimated objective function is

$\hat{J}\left(\theta \right)=\mathrm{ReLU}\left(\underset{z~P_z}{𝔼}\left[\hat{D}\left(G\left(z;\theta \right)\right)\right]-\left(1-\eta \right)\right)$

• • where {circumflex over (D)}(ϕ)=D(τ_E, τ_S, {circumflex over (l)}₁(ϕ), {circumflex over (l)}₂(ϕ)). Consequently, as a direct consequence of Proposition 1 and the universal approximation properties in (37), it is straightforward to show that Ĵ satisfies the Lipschitz inequality

$\begin{array}{cc}❘\hat{J}\left({\theta }_1\right)-\hat{J}\left({\theta }_2\right)❘\le L_D\left(L_1+L_2+\frac{{ϵ}_1+{ϵ}_2}{d_g}\right)L_G{\theta }_1-{\theta }_2& \left(40\right)\end{array}$

• • for all θ₁, θ₂∈Θ_g, where d_gmax_{θ1,θ2∈Θg}∥θ₁−θ₂∥ is the diameter of the compact set Θ_g. Following similar arguments as in the proof of Theorem 3,

$\begin{array}{cc}\begin{array}{c}{{\theta }_g^{\left(i+1\right)}-{\theta }_g^{★}}^2={{\theta }_g^{\left(i\right)}-{\lambda }_i{\hat{h}}^{\left(i\right)}-{\theta }_g^{★}}^2\\ \le {{\theta }_g^{\left(i\right)}-{\theta }_g^{★}}^2-2{\lambda }_i\left[J\left({\theta }_g^{\left(i\right)}\right)-J\left({\theta }_g^{★}\right)\right]+{\lambda }_i^2{{\hat{h}}^{\left(i\right)}}^2\end{array}& \left(41\right)\end{array}$ $\mathrm{Thus},$ $\begin{array}{c}{\hat{J}}_{best}^{\left(i\right)}-J\left({\theta }_g^{★}\right)\le \frac{{{{\theta }_g^{\left(0\right)}-{\theta }_g^{★}}^2+{\sum }_{j=0}^i{\lambda }_j^2{\hat{h}}^{\left(j\right)}}^2}{2{\sum }_{j=0}^i{\lambda }_j}\\ \le \frac{{{\theta }_g^{\left(0\right)}-{\theta }_g^{★}}^2+{L_D^2\left(L_1+L_2+\frac{{ϵ}_1+{ϵ}_2}{d_g}\right)}^2{L}_G^{2}M}{2{\sum }_{j=0}^i{\lambda }_j}\\ \le \frac{{{\theta }_g^{\left(0\right)}-{\theta }_g^{★}}^2+{L_D^2\left(L_1+L_2\right)}^2{L}_G^{2}M}{2{\sum }_{j=0}^i{\lambda }_j}+\frac{L_D^2\left(2L_1+2L_2+\frac{{ϵ}_1+{ϵ}_2}{d_g}\right)\left(\frac{{ϵ}_1+{ϵ}_2}{d_g}\right)L_G^{2}M}{2{\sum }_{j=0}^i{\lambda }_j}\\ \le ϵ\mathrm{for}\mathrm{all}i\ge {\hat{N}}_{ϵ}\end{array}$ $\mathrm{Where}:$ ${\hat{N}}_{ϵ}=\min \left\{i\in \mathbb{Z}❘{\sum }_{j=0}^i{\lambda }_j\ge a+b\right\}$ $\mathrm{With}:$ $a=\frac{{{\theta }_g^{\left(0\right)}-{\theta }_g^{★}}^2+{L_D^2\left(L_1+L_2\right)}^2{L}_G^{2}M}{2ϵ}$ $b=\frac{L_D^2\left(2L_1+2L_2+\frac{{ϵ}_1+{ϵ}_2}{d_g}\right)\left(\frac{{ϵ}_1+{ϵ}_2}{d_g}\right)L_G^{2}M}{2ϵ}$ $\mathrm{From}\left(33\right),N_{ϵ}=\min \left\{i\in \mathbb{Z}❘{\sum }_{j=0}^i{\lambda }_j\ge a\right\}.\mathrm{Thus}$ $\begin{array}{c}{\hat{N}}_{ϵ}=\min \left\{i\in \mathbb{Z}❘{\sum }_{j=0}^i{\lambda }_j-{\sum }_{j=0}^{N_{ϵ}}{\lambda }_j\ge a-{\sum }_{j=0}^{N_{ϵ}}{\lambda }_j+b\right\}\\ =\min \left\{i\in \mathbb{Z}❘{\sum }_{j=N_{ϵ}+1}^i{\lambda }_j\ge a-{\sum }_{j=0}^{N_{ϵ}}{\lambda }_j+b\right\}\end{array}$ $\mathrm{Since}{\sum }_{j=0}^{N_{ϵ}}{\lambda }_j\ge a,\mathrm{it}\mathrm{follows}\mathrm{that}{\hat{N}}_{ϵ}\le \min \left\{i\in \mathbb{Z}❘{\sum }_{j=N_{ϵ}+1}^i{\lambda }_j\ge b\right\}={\overline{N}}_{ϵ}$

Essentially, the above theorem shows that, with the uncertainties on the l₁, l₂ functions, the feasible set Θ*_g(η) can still be reached but with a particular amount of increased convergence time. The increase in the convergence time depends on the choice of step size and the uncertainty bounds.

Example Case Studies

The example implementation included case studies for a power grid system and gas pipeline network. The case studies included attack policies that can be used in implementations of the present disclosure.

The study considered a structure of attack policy as:

$\begin{array}{cc}\pi \left({\theta }_a\right)=\left\{t_0\left({\theta }_a\right),T\left({\theta }_a\right),e\left({\theta }_a,t\right)\right\}& \left(42\right)\end{array}$

It includes the start time of the attack injection t₀∈[t_l, t_u], the duration of the attack injection T∈[T_l, T_u], and the attack profile e(θ_a, t). [t_l, t_u] and [T_l, T_u] are the pre-defined interval of the start time and duration of the attack injection respectively. Examples of such attack policies can be found in FIG. 15, such as ramp attack (RA), sine attack (SA), and pulse attack (PA).

The attack injection time t₀ and the duration time T are determined as

$\begin{array}{cc}{t}_0=t_l+{\theta }_{a1}\left(t_u-t_l\right)& \left(43\right)\end{array}$ $\begin{array}{cc}T=T_l+{\theta }_{a2}\left(t_u-t_l\right)& \left(43\right)\end{array}$

The attack profile functions are determined by two parameters e₁∈[e₁, ē₁], e₂∈[e₂, ē₂], where [e₁, ē₁] and [e₂, ē₂] are also

• • the pre-defined intervals. The value of e₁, e₂ are determined by the output of the generator θ_a:

$\begin{array}{cc}{e}_1={\underset{\_}{e}}_1+{\theta }_{a3}\left({\overline{e}}_1-{\underset{\_}{e}}_1\right)& \left(44\right)\end{array}$ $\begin{array}{cc}{e}_2={\underset{¯}{e}}_2+{\theta }_{a4}\left({\overline{e}}_2-{\underset{¯}{e}}_2\right)& \left(44\right)\end{array}$

Then, for RA policy, the attack profile function is given as:

$\begin{array}{cc}{e}_{RA}\left({\theta }_a,t\right)=\min \left\{\frac{e_2}{e_1}\left(t-t_0\right),e_2\right\}\odot \left(t_0\le t\le t_0+T\right)& \left(45\right)\end{array}$

For SA policy, it is given as

$\begin{array}{cc}{e}_{SA}\left({\theta }_a,t\right)=e_1\sin \left(e_{2}t\right)\odot \left(t_0\le t\le t_0+T\right)& \left(46\right)\end{array}$

For PA policy, it is given as

$\begin{array}{cc}{e}_{PA}\left({\theta }_a,j\right)=e_1\mathrm{sgn}\left(\sin \left(e_2\pi t\right)\right)\odot \left(t_0\le t\le t_0+T\right)& \left(47\right)\end{array}$

Power Grid

The study considered a power grid system comprising n_g generator buses and n_l load buses. An example on IEEE 14-bus system (n_g=5, n_l=14) is used, as shown in FIG. 16. The system is equipped with IoT net meters and frequency meters that measure the net power of buses and generators' frequency. A control center collects the sensor readings through a communication network and performs state estimation, bad data detection, and regulation control of the generator's rotor angles and frequency. However, an attacker could potentially compromise the communication channels or spoof the smart meters to tamper with the sensor readings received by the control center. This can result in a malicious modification of the power grid's state, causing unexpected damages.

The study included an example system's dynamical models and the nominal control setup. Without loss of generality, several assumptions are admitted [34A]:

• • 1. The network voltages stay constant at v=1 p.u.;
  • 2. The angular difference between each bus is small;
  • 3. The conductance is negligible, therefore the system is lossless. The active power flow dynamics at bus i is represented as

$\begin{array}{cc}0=-b_{ij}\sin \left({\theta }_i-{\theta }_j\right)-P_i\left(\theta \right)& \left(48\right)\end{array}$

• • where θ_i denotes i-th bus angle, P_i denotes the active power at the bus i, and b_ij is the associated susceptance of the transmission line between node i and node j. Then, given Ψ(θ) defined as in (2) with ψ(θ_i, θ_j)=−b_ij sin(θ_i−θ_j), the networked power flow dynamics is described by

$\begin{array}{cc}0=\left(A\odot \Psi \left(\theta \right)\right)1-P\left(\theta \right)& \left(49\right)\end{array}$

• • where the network power active power vector can be written as

$P\left(\theta \right)=\left[\begin{array}{c}{P}_g^{n\omega }\left(\theta \right)\\ P_l^{n\omega }\left(\theta \right)\end{array}\right],$

and P_g^nω(θ), P_l^nω(θ) are the active power vector at generator buses and load buses, respectively. Furthermore, the nonlinear swing dynamics of generators is described by [34A]:

$\begin{array}{cc}\left[\begin{array}{cc}I& 0\\ 0& M_g\end{array}\right]\left[\begin{array}{c}\stackrel{.}{\delta }\\ \stackrel{.}{\omega }\end{array}\right]=\left[\begin{array}{cc}0& I\\ 0& -D_g\end{array}\right]\left[\begin{array}{c}\delta \\ \omega \end{array}\right]+\left[\begin{array}{c}0\\ P_g-P_g^{n\omega }\left(\theta \right)\end{array}\right],& \left(50\right)\end{array}$

• • where δ∈^ng is the vector of generator rotor angles, ω∈^ng is the generator frequency vector, M_g is a diagonal matrix of inertial constants for each generator and D_g is a diagonal matrix of damping coefficients and P_g∈^ng contains net power injected at the generator buses. In addition, at load buses, the active power P_l^nω(θ) is equal to the injected power P_l, so that

$\begin{array}{cc}{P}_l^{n\omega }\left(\theta \right)=P_l& \left(51\right)\end{array}$

Since the network power dynamics has already been at its steady state, the supervisory controller in (11) is not required, and the reference active power for generators can be determined as

$\begin{array}{cc}P\left(\theta \right)=\left(A\odot \Psi \left(\theta \right)\right)1& \left(52\right)\end{array}$

By ordering the buses such that the generator nodes appear first, and linearizing the model in (52), the study obtained:

$\begin{array}{cc}\left[\begin{array}{c}{P}_g^{n\omega }\left(\theta \right)\\ P_l^{n\omega }\left(\theta \right)\end{array}\right]=-L\theta & \left(53\right)\end{array}$

• • where the admittance-weighted Laplacian matrix is given by as

$L=\left[\begin{array}{cc}{L}_{\mathrm{gg}}& L_{\lg }\\ L_{\mathrm{gl}}& L_{ll}\end{array}\right]\in {ℝ}^{N\times N},N=n_g+n_l,\mathrm{and}\theta =\left[\begin{array}{c}\delta \\ {\theta }_l\end{array}\right]$

in which θ_l is the bus angle of load buses. Combining (51) and (53) yields

$\begin{array}{cc}{\theta }_l=-L_{ll}^{-1}\left(L_{lg}\delta -P_l\right)& \left(54\right)\end{array}$ $\mathrm{thus},$ $P_g^{n\omega }\left(\theta \right)=-L_{\mathrm{gg}}\delta -L_{lg}{\theta }_l$

Then substituting (55) into the generator dynamical equation (50) yields

$\begin{array}{cc}\stackrel{.}{x}=\underset{A^{\prime }}{\underbrace{\left[\begin{array}{cc}0& I\\ -M_g^{-1}\left(L_{\mathrm{gg}}-L_{\mathrm{gl}}{L}_{\mathrm{ll}}^{-1}{L}_{\lg }\right)& M_g^{-1}{D}_g\end{array}\right]}}x+\underset{B}{\underbrace{\left[\begin{array}{cc}0& 0\\ M_g^{-1}& -M_g^{-1}{L}_{\mathrm{gl}}{L}_{\mathrm{ll}}^{-1}\end{array}\right]}}u& \left(56\right)\end{array}$

• • where the state variables x=[δ^Tω^T]^T∈^2ng contain the generator rotor angles and the generator frequency. The control input vector u=[P_g^T P_l^T]^T∈^ng+nl contains the power injected at the generator buses and load buses. The study assumes the generator frequency and the net power at load buses are measured, then the following measurement model is utilized:

$\begin{array}{cc}y=\underset{C}{\underbrace{\left[\begin{array}{cc}0& I\\ -P\left(\theta \right)L_{\mathrm{ll}}^{-1}{L}_{\lg }& 0\end{array}\right]}}x+\underset{D}{\underbrace{\left[\begin{array}{cc}0& 0\\ P\left(\theta \right)L_{\mathrm{ll}}^{-1}& 0\end{array}\right]}}u& \left(57\right)\end{array}$

• • where y=[ω^T P_net^T]^T∈^ng+nl contain the generator frequency and the net power injected at the load buses. With a sample time T_s=0.01 s, the study obtained the discrete-time model of (56) and (57):

$x\left(k+1\right)=A_{d}x\left(k\right)+B_{d}u\left(k\right)$

• • where A_d=A′T_s−I and B_d=BT_s. A Luenberger observer is employed to estimate the states x(k+1) based on the discrete model in (58):

$\hat{x}\left(k+1\right)=A_d\hat{x}\left(k\right)+B_{d}u\left(k\right)+L\left(y\left(k\right)-C\hat{x}\left(k\right)-Du\left(k\right)\right),$

where L is the observer gain satisfying that A_d−LC has all the eigenvalues inside the unit circle.

Then the PI controller:

$\begin{array}{cc}u\left(k\right)=K_P\tilde{x}\left(k\right)+K_I{\int }_0^k\tilde{x}\left(\tau \right)d\tau & \left(60\right)\end{array}$

• • is used to regulate the generators' rotor angles and frequencies. Namely, select gain matrices K_I and K_p to render the matrix

$\begin{array}{cc}0& I\\ {\mathrm{BK}}_I& A^{\prime }+{\mathrm{BK}}_P\end{array}]$

Hurwitz.

Next, the study discussed the implementation of the proposed attack generator on this IEEE 14-bus system. The attacks are injected through the measurements. The total simulation time is set as 8 s, thereby T_sim=8/T_s=800 time steps in total. The study choose the ratio of the mean control error of the generator rotor angles to the maximum nominal³ generator angles to be the effectiveness index

$\begin{array}{cc}{a}_E=l_1\left(e\right)=\frac{\frac{1}{T_{sim}}{\sum }_{k=0}^{T_{\mathrm{sim}}}{\delta \left(k\right)-{\delta }^r\left(k\right)}_2}{\underset{0\le k\le T_{\mathrm{sim}}}{\max }{{\delta }_{\mathrm{nominal}}\left(k\right)}_2}& \left(61\right)\end{array}$

• • and choose the maximum estimation residual to be the stealthiness index

$\begin{array}{cc}{a}_S=l_2\left(e\right)=\underset{0\le k\le T_{\mathrm{sim}}}{\max }{y\left(k\right)-Du\left(k\right)-C\hat{x}\left(k\right)}_2& \left(62\right)\end{array}$

The three attack policies, shown in FIG. 15, are all tested. The study used the same pre-defined range of start time and duration time of the attack injection t₀∈[0.1 T_sim, 0.2 T_sim], T∈[1e-5, T_sim]. The study choose e₁∈[0.1 T_sim, 0.6 T_sim], e₂∈[−0.04, 0.04] for ramp attack policy, e₁∈[−0.04,0.04], e₂∈[−5,5] for sine attack policy, e₁∈[−0.04,0.04], e₂∈[−5,5] for pulse attack policy. The generator has 4 layers, namely ReLU, ReLU, Tanh, and Sigmoid with 500, 1000, 500 and 4*(n_g+n_l) neurons respectively, respectively. The stealth net includes of 3 ReLU hidden layers and a Sigmoid output layer. The effect net includes of 3 ReLU hidden layers and a Linear output layer. Both the effect net and stealth net have 500, 1000, and 500 neurons at the hidden layers. During each training epoch, they were trained in 20 batches with a size of 2000, and the generator was trained in 10 batches with a size of 5000.

All three attack policies are used separately in the training process of Algorithm 1 shown in FIG. 3B. Algorithm 1 Training generative model with data-driven {circumflex over (l)}₁, {circumflex over (l)}₂ Hyperparameters: τ_E (effectiveness threshold), τ_S (stealthiness threshold), η (success probability), π (attack policy), s (scale factor of leakyReLU function), P_z (prior distribution, i.e. uniform distribution), N_sample (number of sample points for each epoch), N_epoch (total number of epochs for the training)

The same stealthiness and effectiveness thresholds are chosen: ϵ=0.05 and α=49%, so the feasible attack set is defined as

$\begin{array}{cc}=\left(0.05,49\%\right)=\left\{e❘l_1\left(e\right)\ge 49\%,l_2\left(e\right)\le 0.05\right\}& \left(63\right)\end{array}$

For Algorithm 1, other hyperparameters are chosen as η=0.8 and s=0. FIG. 17, FIG. 18, and FIG. 19 depict the testing results of the trained attack generator in the system simulation after each epoch for the three attack policies. These results indicate the successful training of the attack generation framework. Initially, most of the generated attacks could be detected and could not arouse effect. So the generator tried to decrease the stealthiness, but doing so also affected the effectiveness because of the positive correlation between the two metrics. Once the stealthiness went below the threshold, the generator began exploring more effective attacks until the effectiveness exceeded the given threshold. By comparing the testing performance policy-wise, it is seen that the attack generation framework approached the boundary of stealthiness more and achieved bigger effectiveness with the SA and PA policies compared to the RA policy. This finding implies that SA and PA policies enable more vulnerability space, which the attack generation framework could explore, than the RA policy. Additionally, FIG. 29 presents the corresponding testing scores, which indicate the ratio of feasible attacks in (0.05, 49%) to the total testing samples for the epochs. It shows that the generator with SA policy converges faster than the other two, and the generator with PA policy is the slowest, likely due to its complexity.

Next, the study evaluated the performance of the generated attacks with different attack policies by injecting them into the system's time-series simulation. As shown in FIG. 29, the study used the trained attack generators at epoch 10 for RA policy and epoch 8 for SA and PA policies. The examples of the time-series performance of the generated attacks are shown in FIG. 20. The first row displays the attack signals generated by the trained attack generators for all 19 sensors. The second to sixth rows show the comparison of generator rotor angles before and after the attack injection. The last row shows the outputted estimation residual (stealthiness index) by BDD. The results show that the proposed attack generation framework with PA policy could explore the most vulnerable space of the system since the largest BDD outputs touch the detection boundary during the simulation. Furthermore, different attack policies produce different effects on the system. RA policy introduces a constant bias on sensor readings, and thus the shapes of the attacked rotor angles are similar to the nominal ones. SA policy introduces waves in the rotor angles, while PA policy introduces higher-frequency fluctuations in rotor angles. Nonetheless, the proposed attack generation framework can find the corresponding feasible attack set for any attack policy.

Gas Pipeline

The study further considered gas pipeline systems equipped with IoT sensors and actuators, a supervisory controller, and automatic controllers, as shown in FIG. 21.

A nonlinear creep flow model is used to describe the gas flow in a pipeline segment [35A]:

$\frac{\partial p}{\partial t}=-\frac{c^2\rho }{G}·\frac{\partial Q_n}{\partial x}$

• • where p is the pressure in the pipeline, x is the position of flow, c is the speed of sound in gas [m/s], ρ is the average gas density over the cross-section area of the pipeline [kg/m³], Q_n is the volumetric flow at standard conditions, f is the friction factor, D is the diameter of the pipeline and G=π(D/2)². While the partial differential equation (PDE) model provides an accurate description of the gas flow dynamics for infrequent online optimization, an ordinary differential equation (ODE) model is sufficient to describe the pressure dynamics at nodes [35A]. The dynamical pressure model at node i is given as

$\begin{array}{cc}{\dot{p}}_i=\frac{c^2}{V_i}\sum _{j\in {𝒩}_i}\sqrt{\frac{❘p_j^2-p_i^2❘}{K_{ij}}}\mathrm{sgn}\left(p_j-p_i\right)-w_i& \left(65\right)\end{array}$

• • where p_i is the pressure at node i, N_i denote the set of neighborhood pipelines of node i,

$K_{\mathrm{ij}}=\frac{64f_{\mathrm{ij}}{c}^2\Delta x_{\mathrm{ij}}}{{\pi }^2{D}_{\mathrm{ij}}^5},$ $V_i=\frac{\pi }{8}\sum _{j\in N_i}{D}_{\mathrm{ij}}^2\Delta x_{\mathrm{ij}}$

• • and w_i is the mass flow at node i, Δx_ij denotes the length of pipeline between node i and node j. Here,

$\begin{array}{cc}\psi \left(p_i,p_j\right)=\sqrt{\frac{❘p_j^2-p_i^2❘}{K_{\mathrm{ij}}}}\mathrm{sgn}\left(p_j-p_i\right)& \left(66\right)\end{array}$

• • then the model of the entire pipeline network is given by

$\stackrel{.}{p}=c^2{V}^{-1}\left(A\odot \Psi \left(p\right)\right)1-w$

• • where Ψ(θ) is defined as in (2) with ψ(θ_i, θ_j) defined in (66), V⁻¹=diag([V₁⁻¹, V₂⁻¹, . . . , V_n⁻¹]), w_sup is a vector of mass inflows at the supply points and w_dem is a vector of mass outflows at the demand points. Clearly w=−B_supw_sup+B_demw_dem, where B_dem and B_sup are appropriately dimensioned matrices such that B_sup^TB_dem=0 and P[−B_sup|B_dem]=I for some projection matrix P.

In addition, a dynamical model of the compressor between two nodes i and j is given by [36A]:

${\dot{\omega }}_{\mathrm{ij}}=\frac{1}{J_{o_i}}\left(u_{\mathrm{ij}}+\rho r_2^2{q}_{\mathrm{ij}}{\omega }_{\mathrm{ij}}\right)$ ${\dot{w}}_{\mathrm{ij}}=q_{\mathrm{ij}}$ $\mathrm{where}$ $\varphi \left(\omega ,q\right)={\left(1+\frac{1}{T_{\mathrm{in}}{c}_{p}1000}\left(\rho r_2^2{\omega }^2+k_f{q}^2-\frac{r_1^2}{2}{\left(\omega -\frac{q}{A_c{r}_1{\rho }_{\mathrm{in}}}\right)}^2\right)\right)}^{\frac{\gamma }{\gamma -1}},$

ω_ij is the compressor rotor angular velocity [rad/s], q_ij is the mass power flow [kg/s], ϕ(ω, q) is the compressor characterization, u_ij is the mechanical torque [N-m] applied to the rotor shaft or inertia J_0i [kg-m²]. p_i and p_j are the input and output pressures of the compressor, respectively. For the system considered, the compressors are at the extremes of the network. Thus, at each compressor node, one of p_i and p_j will be the ambient pressure. p_i=P_amb for upstream compressors, while p_j=P_amb for downstream compressors. The compressor could be controlled by a simple PID controller:

$\begin{array}{cc}{u}_{\mathrm{ij}}=K_P{\tilde{w}}_{\mathrm{ij}}\left(t\right)+K_I{\int }_0^t{\tilde{w}}_{\mathrm{ij}}\left(\tau \right)d\tau +K_D\frac{d{\tilde{w}}_{\mathrm{ij}}\left(t\right)}{\mathrm{dt}}& \left(69\right)\end{array}$

• • where {tilde over (w)}_ij(t)=w_ij−w_ij^r. Similarly, the PID gains K_P, K_I, K_D are designed to achieve the closed-loop performance in (7). Linearizing the models (67) around the equilibrium point (p_eq, w_eq) and discretizing with a fixed time step T_s=0.01 s yields

$\begin{array}{cc}{\tilde{p}}_{k+1}=A_d{\tilde{p}}_k+B_d{\tilde{w}}_k& \left(70\right)\end{array}$

• • where {tilde over (p)}_k=p_k−p_eq, {tilde over (w)}_k=w_k−w_eq, and A_d, B_d are given as in (10) with

$$\begin{gathered} \frac{\partial \psi \left(p_i,p_j\right)}{{\partial }_{p_i}}=-\frac{p_i}{K_{\mathrm{ij}}}\sqrt{\frac{K_{\mathrm{ij}}}{❘p_j^2-p_i^2|}}\mathrm{sgn}\left(p_j-p_i\right)\mathrm{sgn}\left(p_j^2-p_i^2\right),\text{ \\ ∂ψ⁡(pi,pj)∂pj=piKij⁢Kij❘"[LeftBracketingBar]"pj2-pi2|⁢sgn⁢(pj-pi)⁢sgn⁢(pj2-pi2).} \end{gathered}$$

The supervisory target generator is given in (11) to regulate the node pressures around its equilibrium point at a lower frequency. The study assumed all node pressures and mass flows could be directly measured subject to noise v:

$\begin{array}{cc}y=x^{\prime }+v& \left(71\right)\end{array}$ $\mathrm{where}$ $x^{\prime }\stackrel{△}{=}\left[\begin{array}{c}p\\ w_{\sup }\\ w_{\mathrm{dem}}\end{array}\right]$

composes of the node pressure states, the wellhead supplied mass flow and the demand mass flow. An unscented Kalman filter (UKF) [37A] is used to estimate the states from the noisy measurements. The study tested three different pipeline systems with distinct network topologies, namely the Linear topology, Tree topology, and Cyclic topology, as illustrated in FIG. 22 [38A]. The linear topology example includes 9 nodes, the tree topology has 11 nodes and the cyclic topology has 13 nodes. All three pipeline systems have two gas supply wellheads, one demand node, and one supply node that requires control.

Next, the study discusses the implementation of the proposed attack generation system in these pipeline systems. The attacks are injected through the measurements, then the measurement model in (71) becomes

$\begin{array}{cc}y=x^{\prime }+v+e& \left(72\right)\end{array}$

The total simulation time is set as 200 s, thereby T_sim=200/T_s=20000 time steps in total. To train the attack generator, the example implementation used the ratio of the mean value of the node pressure control error to a constant nominal pressure as the effectiveness metric:

$\begin{array}{cc}{a}_E=l_1\left(e\right)=\frac{1}{T_{\mathrm{sim}}}\sum _{k=0}^{T_{\mathrm{sim}}}{\tilde{p}\left(k\right)}_2& \left(73\right)\end{array}$

• • and use the maximal estimation residual as the stealthiness metric:

$\begin{array}{cc}{a}_S=l_2\left(e\right)=\underset{0\le k\le T_{\mathrm{sim}}}{\max }{y\left(k\right)-x^{\prime }\left(k\right)}_2& \left(74\right)\end{array}$

This case study focused on testing the Ramp attack policy, as defined in (45), with the pre-defined range of injection start time t₀∈[0.1 T_sim, 0.2 T_sim], injection duration time T∈[1e-5, 0.5 T_sim] and two attack value function parameters e₁∈[0.1 T_sim, T_sim] and e₂∈[0,0.5*y]. The input of the attack value function (45) is given in (43), where the input of the attack policy g is the output of the generator. For linear topological pipeline system, it has 12 measurements, so the output size of the associated generator is 4*12=48. Similarly, the output size of the generator for the tree topological pipeline system is 4*14=56 and the output size of the generator for the cyclic topological pipeline system is 4*16=64. The input sizes of the generators are the same and equal to 10, and they all have 4 layers: ReLU (500 neurons), Re LU (1000 neurons), Tanh (500 neurons), Sigmoid. The stealth net includes of 3 ReLU hidden layers and a Sigmoid output layer, and the effect net includes of 3 ReLU hidden layers and a Linear output layer. The numbers of neurons at hidden layers are 500, 1000, 500 for both the effect net and stealth net. In each training epoch, the discriminators are trained in 20 batches with a size of 1000, and the generator is trained in 5 batches with a size of 5000. The same stealthiness and effectiveness thresholds are chosen for all three topological pipeline systems: ϵ=0.02 and α=65, so the feasible attack set is defined as:

$\begin{array}{cc}=\left(0.05,65\%\right)=\left\{e❘l_1\left(e\right)\ge 65\%,l_2\left(e\right)\le 0.05\right\}& \left(75\right)\end{array}$

For Algorithm 1, other hyperparameters are chosen as η=0.8 and s=0. All generators converge in the first 2 epochs. Therefore, in order to show the convergence process of the generative system, the study plotted the testing results after each batch of training in the first 2 epochs.

FIG. 23, FIG. 24, and FIG. 25 show the testing performance of the 180 sample generative attacks from the trained generator after each batch in the first 2 epochs, and FIG. 30 shows the testing scores of the trained generators after batches in the first 2 epochs. As shown in FIG. 23, FIG. 24, and FIG. 25, it is seen that all three generators started at a similar situation where although the generative attacks are not detectable, they are not effective enough. Thus, the generators attempted to increase the effectiveness in the first epoch, while stealthiness also increased due to the positive correlation between effectiveness and stealthiness. For the linear topological case as shown in FIG. 23, the generator converged at the second epoch since the effectiveness and stealthiness do not move anymore. However, a small portion of generative attacks is detectable. This is because the threshold setting is a little strict for this case. If a lower effectiveness threshold is given, the testing score of the generator becomes 100% after converging. For the tree topological case as shown in FIG. 24, the generator was also too aggressive to increase effectiveness, which caused the stealthiness to exceed the thresholds in the first epoch. But in the second epoch, the generator was able to lower it back and converged to the target feasible attack set. For the cyclic topological case as shown in FIG. 25, the generator could converge to the biggest effectiveness while maintaining the smallest stealthiness compared other two topological cases. This indicates the cyclic topological pipeline network used in this case study is the most vulnerable one and the linear topological pipeline network used in this case study is the most resilient one.

Next, the study evaluated the performance of generated feasible attacks by injecting them into the time-series simulation of those three pipeline systems. The results are presented in FIG. 26, FIG. 27, and FIG. 28. The subplots show the ratio of the attack signals e to the measurements. The attacks use Ramp attack policy injected between 0.1 T_sim=20 s and 0.2 T_sim=40 s. The subplots show the BDD output which is the estimation residual ∥y(k)−x′(k)∥₂. The subplots further show the node pressure deviation ratio from the equilibrium points. The generative attacks misled the node pressure to deviate from the target equilibrium point and arrive at a wrong equilibrium point. Comparing the performance topology-wise, it is seen that the most deviation from the nominal node pressure is aroused for the cyclic topological network. This phenomenon is consistent with the testing results shown in FIG. 23, FIG. 24, and FIG. 25. It indicates that more nodes and a more complicated network topology would provide more opportunities for the proposed attack generation framework to explore more vulnerability space.

The study developed a data-driven attack generation system accessing the vulnerability of nonlinear NCPS. It does not require the model knowledge of the system but only runtime data. The framework comprises three interactive models: two discriminative models that assess attack signals and a generative model that generates feasible attacks. To ensure a high probability of successful attack generation, the study disclosed a loss function and a training algorithm with convergence guarantee. The implementations described herein are highly versatile and can be applied to various NCPS. It can enable the learning of stealthiness boundaries while exploring effectiveness, which holds immense potential for use in resilient co-design with learning-based BDDs.

The study further contemplates that implementations of the present disclosure can automatically find an optimal attack policy using an algorithm similar to policy learning and/or include an automatic approach to identify the thresholds of effectiveness and stealthiness, ensuring that the feasible attack set is not void. Implementations of the present disclosure can include additional features not described as part of algorithm 1 shown in FIG. 3B, for example. Additionally, it should again be understood that the example applications described in the study are non-limiting examples, and that any number/type of systems can be studied using the disclosed attack generation systems.

REFERENCES

Although the subject matter has been described in language specific to structural features and/or methodological acts, it is to be understood that the subject matter defined in the appended claims is not necessarily limited to the specific features or acts described above. Rather, the specific features and acts described above are disclosed as example forms of implementing the claims.

The following patents, applications, and publications, as listed below and throughout this document, describes various application and systems that could be used in combination the exemplary system and are hereby incorporated by reference in their entirety herein

• • [1] A. Hobbs, “The colonial pipeline hack: Exposing vulnerabilities in us cybersecurity,” in SAGE Business Cases. SAGE Publications: SAGE Business Cases Originals, 2021.
  • [2] R. Langner, “Stuxnet: Dissecting a cyberwarfare weapon,” IEEE Security & Privacy, vol. 9, no. 3, pp. 49-51, 2011.
  • [3] F. Pasqualetti, F. Dorfler, and F. Bullo, “Attack detection and identification in cyber-physical systems,” IEEE transactions on automatic control, vol. 58, no. 11, pp. 2715-2729, 2013.
  • [4] G. Chen, Y. Zhang, S. Gu, and W. Hu, “Resilient state estimation and control of cyber-physical systems against false data injection attacks on both actuator and sensors,” IEEE Transactions on Control of Network Systems, vol. 9, no. 1, pp. 500-510, 2021.
  • [5] T. Sui, Y. Mo, D. Marelli, X. Sun, and M. Fu, “The vulnerability of cyber-physical system under stealthy attacks,” IEEE Transactions on Automatic Control, vol. 66, no. 2, pp. 637-650, 2020.
  • [6] X. Ren and Y. Mo, “Secure detection: Performance metric and sensor deployment strategy,” IEEE Transactions on Signal Processing, vol. 66, no. 17, pp. 4450-4460, 2018.
  • [7] M. H. Shahriar, A. A. Khalil, M. A. Rahman, M. H. Manshaei, and D. Chen, “iattackgen: Generative synthesis of false data injection attacks in cyber-physical systems,” in 2021 IEEE Conference on Communications and Network Security (CNS). IEEE, 2021, pp. 200208.
  • [8] A. Khazraei, S. Hallyburton, Q. Gao, Y. Wang, and M. Pajic, “Learning-based vulnerability analysis of cyber-physical systems,” in 2022 ACM/IEEE 13th International Conference on Cyber-Physical Systems (ICCPS). IEEE, 2022, pp. 259-269.
  • [9] A. N. Kolmogorov and A. T. Bharucha-Reid, Foundations of the theory of probability: Second English Edition. Courier Dover Publications, 2018.
  • [10] C. Godsil and G. F. Royle, Algebraic graph theory. Springer Science & Business Media, 2001, vol. 207.
  • [11] A. Osiadacz, Simulation and analysis of gas networks. Gulf Publishing Company, Houston, TX, 1987.
  • [12] H. Perez-Blanco and T. B. Henricks, “A gas turbine dynamic model for simulation and control,” in Turbo Expo: Power for Land, Sea, and Air, vol. 78637. American Society of Mechanical Engineers, 1998, p. V002T03A002.
  • [13] Y. Liu, P. Ning, and M. K. Reiter, “False data injection attacks against state estimation in electric power grids,” ACM Transactions on Information and System Security (TISSEC), vol. 14, no. 1, pp. 133, 2011.
  • [14] Y. Zheng and O. M. Anubi, “Resilient observer design for cyberphysical systems with data-driven measurement pruning,” in Security and Resilience in Cyber-Physical Systems. Switzerland: Springer, 2022, pp. 85-117.
  • [15] Y. Mo and B. Sinopoli, “False data injection attacks in control systems,” in Preprints of the 1st workshop on Secure Control Systems, vol. 1, 2010.
  • [16] S. Zhao, H. Ren, A. Yuan, J. Song, N. Goodman, and S. Ermon, “Bias and generalization in deep generative models: An empirical study,” Advances in Neural Information Processing Systems, vol. 31, 2018.
  • [1A] M. Dahleh, E. Frazzoli, A. Megretski, S. Ozdaglar, P. Parrilo, and D. Shah, “Information and control: a bridge between the physical and decision layers,” in NSF Workshop on CyberPhysical Systems, 2006.
  • [2A] Q. Zhu and L. Bushnell, “Networked cyber-physical systems: Interdependence, resilience and information exchange,” in 2013 51st Annual Allerton conference on communication, control, and computing (Allerton). IEEE, 2013, pp. 763-769.
  • [3A] S. Paul, Z. Ni, and C. Mu, “A learning-based solution for an adversarial repeated game in cyber-physical power systems,” IEEE Transactions on Neural Networks and Learning Systems, vol. 31, no. 11, pp. 4512-4523, 2019.
  • [4A] C. Spandonidis, P. Theodoropoulos, F. Giannopoulos, N. Galiatsatos, and A. Petsa, “Evaluation of deep learning approaches for oil & gas pipeline leak detection using wireless sensor networks,” Engineering Applications of Artificial Intelligence, vol. 113, p. 104890, 2022.
  • [5A] L. Zhang, J. Sun, and G. Orosz, “Hierarchical design of connected cruise control in the presence of information delays and uncertain vehicle dynamics,” IEEE Transactions on Control Systems Technology, vol. 26, no. 1, pp. 139-150, 2017.
  • [6A] P. P. Sahoo, A. Kanellopoulos, and K. G. Vamvoudakis, “Intermittent learning through operant conditioning for cyber-physical systems,” IEEE Transactions on Neural Networks and Learning Systems, 2021.
  • [7A] F. Pasqualetti, F. Dorfler, and F. Bullo, “Attack detection and identification in cyber-physical systems,” IEEE transactions on automatic control, vol. 58, no. 11, pp. 2715-2729, 2013.
  • [8A] K.-D. Lu, L. Zhou, and Z.-G. Wu, “Representation-learning-based cnn for intelligent attack localization and recovery of cyber-physical power systems,” IEEE Transactions on Neural Networks and Learning Systems, 2023.
  • [9A] Y. Liu, P. Ning, and M. K. Reiter, “False data injection attacks against state estimation in electric power grids,” ACM Transactions on Information and System Security (TISSEC), vol. 14, no. 1, pp. 1-33, 2011.
  • [10A] Y. Mo and B. Sinopoli, “False data injection attacks in control systems,” in Preprints of the 1st workshop on Secure Control Systems, vol. 1, 2010.
  • [11A] L. Hu, Z. Wang, Q.-L. Han, and X. Liu, “State estimation under false data injection attacks: Security analysis and system protection,” Automatica, vol. 87, pp. 176-183, 2018.
  • [12A] T. Sui, Y. Mo, D. Marelli, X. Sun, and M. Fu, “The vulnerability of cyber-physical system under stealthy attacks,” IEEE Transactions on Automatic Control, vol. 66, no. 2, pp. 637-650, 2020.
  • [13A] X. Liu, Z. Bao, D. Lu, and Z. Li, “Modeling of local false data injection attacks with reduced network information,” IEEE Transactions on Smart Grid, vol. 6, no. 4, pp. 1686-1696, 2015.
  • [14A] A.-Y. Lu and G.-H. Yang, “False data injection attacks against state estimation without knowledge of estimators,” IEEE Transactions on Automatic Control, 2022.
  • [15A] A. Khazraei, S. Hallyburton, Q. Gao, Y. Wang, and M. Pajic, “Learning-based vulnerability analysis of cyber-physical systems,” in 2022 ACM/IEEE 13th International Conference on Cyber-Physical Systems (ICCPS). IEEE, 2022, pp. 259-269.
  • [16A] I. L. Carreno, A. Scaglione, A. Zlotnik, D. Deka, and K. Sundar, “An adversarial model for attack vector vulnerability analysis on power and gas delivery operations,” Electric Power Systems Research, vol. 189, p. 106777, 2020.
  • [17A] R. Fu, X. Huang, Y. Xue, Y. Wu, Y. Tang, and D. Yue, “Security assessment for cyber physical distribution power system under intrusion attacks,” IEEE Access, vol. 7, pp. 75615-75628,2018.
  • [18A] K. Yan, X. Liu, Y. Lu, and F. Qin, “A cyber-physical power system risk assessment model against cyberattacks,” IEEE Systems Journal, 2022.
  • [19A] M. Mohammadpourfard, F. Ghanaatpishe, M. Mohammadi, S. Lakshminarayana, and M. Pechenizkiy, “Generation of false data injection attacks using conditional generative adversarial networks,” in 2020 IEEE PES Innovative Smart Grid Technologies Europe (ISGT-Europe). IEEE, 2020, pp. 41-45.
  • [20A] M. H. Shahriar, A. A. Khalil, M. A. Rahman, M. H. Manshaei, and D. Chen, “iattackgen: Generative synthesis of false data injection attacks in cyber-physical systems,” in 2021 IEEE Conference on Communications and Network Security (CNS). IEEE, 2021, pp. 200-208.
  • [21A] C. Godsil and G. F. Royle, Algebraic graph theory. Springer Science & Business Media, 2001, vol. 207.
  • [22A] G. F. Coulouris, J. Dollimore, and T. Kindberg, Distributed systems: concepts and design. pearson education, 2005.
  • [23A] J. K. Hedrick and A. Girard, “Control of nonlinear dynamic systems: Theory and applications,” Controllability and observability of Nonlinear Systems, vol. 48, pp. 133-160, 2005.
  • [24A] M. Pajic, J. Weimer, N. Bezzo, O. Sokolsky, G. J. Pappas, and I. Lee, “Design and implementation of attack-resilient cyberphysical systems: With a focus on attack-resilient state estimators,” IEEE Control Systems Magazine, vol. 37, no. 2, pp. 66-81, 2017.
  • [25A] E. D. Sontag, “A ‘universal’ construction of artstein's theorem on nonlinear stabilization,” Systems & control letters, vol. 13, no. 2, pp. 117-123, 1989.
  • [26A]—, Mathematical control theory: deterministic finite dimensional systems. Springer Science & Business Media, 2013, vol. 6.
  • [27A] B. W. Bequette, Process control: modeling, design, and simulation. Prentice Hall Professional, 2003.
  • [28A] Y. Zheng and O. M. Anubi, “Attack-resilient weighted l₁ observer with prior pruning,” in American Control Conference, 2021, pp. 340-345.
  • [29A] Y. Zheng, S. B. Mudhangulla, and O. M. Anubi, “Moving-horizon false data injection attack design against cyber-physical systems,” Control Engineering Practice, vol. 136, p. 105552, 2023.
  • [30A] C. Huang, J. Wen, Y. Xu, Q. Jiang, J. Yang, Y. Wang, and D. Zhang, “Self-supervised attentive generative adversarial networks for video anomaly detection,” IEEE transactions on neural networks and learning systems, 2022
  • [31A] K. Hornik, M. Stinchcombe, and H. White, “Multilayer feedforward networks are universal approximators,” Neural networks, vol. 2, no. 5, pp. 359-366,1989.
  • [32A] I. Gühring, G. Kutyniok, and P. Petersen, “Error bounds for approximations with deep relu neural networks in w s, p norms,” Analysis and Applications, vol. 18, no. 05, pp. 803-859, 2020.
  • [33A] K. Hornik, M. Stinchcombe, and H. White, “Universal approximation of an unknown mapping and its derivatives using multilayer feedforward networks,” Neural networks, vol. 3, no. 5, pp. 551-560, 1990.
  • [34A] E. Scholtz, “Observer-based monitors and distributed wave controllers for electromechanical disturbances in power systems,” Ph.D. dissertation, Massachusetts Institute of Technology, 2004.
  • [35A] A. Osiadacz, Simulation and analysis of gas networks. Gulf Publishing Company, Houston, TX, 1987.
  • [36A] H. Perez-Blanco and T. B. Henricks, “A gas turbine dynamic model for simulation and control,” in Turbo Expo: Power for Land, Sea, and Air, vol. 78637. American Society of Mechanical Engineers, 1998, p. V002T03A002.
  • [37A] Y. Zheng and O. M. Anubi, “Data-driven vulnerability analysis of networked pipeline system,” in 2023 IEEE Conference on Control Technology and Applications (CCTA). IEEE, 2023, pp. 1022-1027.
  • [38A] R. Z. Ríos-Mercado and C. Borraz-Sánchez, “Optimization problems in natural gas transportation systems: A state-of-the-art review,” Applied Energy, vol. 147, pp. 536-555, 2015.

Claims

1. A computer-implemented method of training a machine learning model, the method comprising: generating a random parameter vector;generating a random attack dataset;inputting a plurality of samples of the random attack dataset into a generator;generating, by the generator, a generated attack dataset;training a discriminator using the random attack dataset and the generated attack dataset; andtraining the generator using the trained discriminators and a loss function.

2. The computer-implemented method of claim 1, wherein the generated attack dataset comprises an attack policy.

3. The computer-implemented method of claim 2, wherein the attack policy comprises a ramp attack.

4. The computer-implemented method of claim 2, wherein the attack policy comprises a sensor attack.

5. The computer-implemented method of claim 2, wherein the generator comprises a deep neural network.

6. The computer-implemented method of claim 5, wherein the generator is trained with the loss function:

7. A computer-implemented method for training a machine learning model, the method comprising: receiving a generative model parameter vector;generating sample points from a prior distribution of the generative model parameter vector;generating attack policy parameterssimulating simulated attack policy parameters by a system simulationtraining a generative model by the simulated attack policy parameters.

8. The computer-implemented method of claim 7, wherein training the generative model comprises a deep neural network.

9. The computer-implemented method of claim 7, wherein training the generative model comprises selecting a best loss value by Jbest(i+1)=min{Jbest(i),J(θg(i+1))}.

10. The computer-implemented method of claim 7, wherein the attack policy parameters comprise a ramp attack policy.

11. The computer-implemented method of claim 7, wherein the attack policy parameters comprise a sine attack policy.

12. The computer-implemented method of claim 7, wherein the attack policy parameters comprise a pulse attack policy.

13. A cybersecurity system comprising: a physical planta network configured for communications and/or control of the physical plant;a cybersecurity controller operably connected to the network, wherein the cybersecurity controller comprises a processor and memory with instructions stored thereon, that, when executed by the processor cause the processor to: receive a trained attack generative model;simulate, by the trained attack generative model, a plurality of attacks on the physical plant;determine an attack effectiveness of each of the plurality of attacks; andcontrol access to the network based on the effectiveness of each of the plurality of attacks.

14. The system of claim 13, wherein the physical plant comprises a networked pipeline system.

15. The system of claim 14, wherein the networked pipeline system comprises a plurality of pressure sensors operably coupled to the network, and controlling access to the network comprises securing at least one of the plurality of pressure sensors.

16. The system of claim 13, wherein the physical plant comprises a power grid.

17. The system of claim 16, wherein the power grid comprises a plurality of meters configured to measure the power and frequency of the power grid, and controlling access to the network comprises securing at least one of the plurality of meters.

18. The system of claim 13, wherein the physical plant comprises a cyber-physical system.

19. The system of claim 13, wherein the physical plant comprises an industrial facility.

20. The system of claim 13, wherein the controller further contains instructions that cause the processor to simulate a vulnerability of the physical plant and control the physical plant based on the vulnerability.