# Toolgraph Design of Optimal and Feasible Control Strategies for Time-Varying Dynamical Systems

Authors: Z. Kowalczuk; K. E. Olinski
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

The paper presents a new method for designing optimal and feasible control strategies for time-variant dynamical processes. The key point of the presented idea lies in utilizing a flow graph structure for representing pertinent properties of the autonomous dynamics of a given dynamical process in a time-and-state space, which is composed of certain elementary segments. The structure is referred to as a time-and-state space toolgraph. In the procedure, each segment of the temporary state space is assigned a node of the time-and-state space toolgraph. The flow values are proportional to the cost of driving the operational point of the dynamical process between the centers of adjacent segments. Any of the discrete optimization algorithms can be applied to search for a cheapest path connecting the initial and terminal points of the sought optimal piecewise-linear trajectory of the operational points in the considered time-and-state space. Additional assumptions or restrictions concerning arbitrary forbidden zones for the operational points can be easily taken into account. In such cases the nodes representing the segments partially or entirely belonging to the finite forbidden zones are deposed from the toolgraph structure.


---

## Body

## 1. Introduction

The general idea of representing a continuous workspace by means of a graph structure has recently been explored, for instance, for searching optimal trajectories of robot movements in 2D [1] and 3D [2] workspaces, as well as for designing optimal control in multidimensional spaces [3].The presented method is a consequent extension of our research track of designing optimal and feasible control by means of discrete optimization methods [4–6]. In particular, certain deterministic time variations of the controlled process are taken into consideration by extending the common state space with an additional time dimension.The algorithm consists in finding a trajectory in such an extended state space (called henceforth the time-state space), which should satisfy two principal prerequisites. First, the sought trajectory circumvents certain finite areas (called the forbidden zones) which can, for instance, represent some failures previously identified. Second, the autonomous dynamics of the process is used in the process of minimizing the adopted indicator of control cost, which is a part of the problem definition and an input to a corresponding minimization procedure. The procedure starts with segmenting the defined operational workspace, which is a subset of the time-and-state space taken into consideration during the process of optimal trajectory design. For each segment a set of representative values is determined. The representative values (representatives) are all the quantities which are necessary to calculate the cost of an operational point movement in a given segment area. The next step is based on a formulation of a flow graph structure (called toolgraph), which reflects the arrangement of the segments along with their appropriate representatives. Any discrete optimization algorithm that performs optimal search for the cheapest path can be utilized for searching the time-and-state space toolgraph for the cheapest path connecting the initial and terminal nodes representing the segments, which contain initial and terminal points of the sought trajectory, respectively.This paper is organized as follows: Section2 gives a formal description of the presented idea and in Section 3 an extended example is presented. The paper ends by the Conclusions section.


## 2. Problem and Solution Description

This section provides some basic definitions (confer [5]) and gives a raw description of the presented method.


### 2.1. Problem Statement

Definition 2.1 (autonomous and forced dynamics).
A dynamical and time-variant process is described by the following function:(2.1)x˙(t)=f(x(t),u(t),t),(2.2)f:ℛl×ℛm×ℛ→ℛl,
where x=[x1,…,xl] is a state vector and u=[u1,…,um] is a control vector, and l,m∈𝒩. Let the autonomous dynamics be described by
(2.3)x˙x(t)=fx(x(t),t)
and let the forced dynamics be defined by a simple centroaffine transformation
(2.4)x˙u(t)=fu(x(t),u(t))=Fu(t),
where F is a fixed (time-invariant) l×m input control matrix, and the above contributors satisfy the following additivity condition:
(2.5)x˙(t)=x˙x(t)+x˙u(t),
where xx and xu are two, autonomous and forced state components.Definition 2.2 (operational subspaceP).
A subsetP=D1×D2×⋯×Dl, where Di, for i=1,…,l, represents a range along the axis xi, which is taken into account while seeking an optimal trajectory, is said to be an operational subspace.Definition 2.3 (forbidden zoneZ).
A subsetZ of P prohibited for operational points is referred to as a forbidden zone. This means that the sought optimal trajectory cannot enter it.Definition 2.4 (system trajectoryE(u)).
Any sequenceE(u) of consecutive states of a given dynamical system for a feasible control input u from a set 𝒰 is said to be a system’s trajectory or a trajectory of operational points in the system state space P∖Z.Definition 2.5 (transition vectorΛ).
A transition vectorΛ is an ordered set of two elements {x0,xk}(2.6)Λ={(x0,xk):x0,xk∈E(u)},
where x0 is the first element and xk is the last element of the sought (optimal) trajectory.Definition 2.6 (cost functionJ(E(u)) of the system and control trajectory).
Let𝒜 be a set of all trajectories E∈𝒜⊂P∖Z. Any real function of the class 𝒜→ℛ+ is said to be a cost function J(E(u)) of this trajectory.Definition 2.7 (optimal trajectoryE*).
LetΞ⊂𝒜 be the subset of all possible trajectories that perform the desired transition Λ (i.e., which start at x0 and terminate at xk). We say that a trajectory E* is (costly) optimal if it satisfies the following conditions:
(2.7)∀E∈ΞJ(E)≥J(E*),∀x∈E*x∈P∖Z.Definition 2.8 (autonomous dynamics mapM).
The imagex˙x of the whole space P in the transformation fx is said to be an autonomous dynamics map M(=x˙x).Our ultimate objective is to find the optimal and feasible controlu*(t) along with its corresponding state trajectory E*=E(u*(t))∈(P∖Z), which implements the complete transition Λ of the dynamical system (2.1) from its initial state x(t0)=x0 to the target state x(tk)=xk and at the same time minimizes the cost functional of the following form:
(2.8)J(E(u))=∫t0tk‍(∑i=1m‍βi|ui(t)|)dt,
where βi, i=1,…,m are nonnegative factors often used for weighting, and tk is a suitable transition-time interval resulting from the applied optimal control procedure.


### 2.2. Solution Method

For the dynamical systems described in Definition2.1, we consider its simple discrete-time variability of the computationally-convenient jump type allowing for the state-space system approximation scheme depicted in Figure 1.Figure 1
Continuous-time dynamic system with discrete-time variability.Leaving apart the initial conditions defined for a certaint0, let us consider consecutive continuous-time computational time moments tn, for n=1,…,N, with tn-tn-1=h=Δt, being a computation step (or a simulation insight), where N∈𝒩 is a discrete-time horizon under consideration.Definition 2.9 (temporary workspacePn).
At each of the above-mentioned computational moments, a separate temporary duplicatePn of the operational subspace P will be associated, as shown in Figure 1.Definition 2.10 (snapshotsMn of the dynamics map).
The multiplication of the subspaceP allows us to define separate snapshots Mn of the dynamics map M assigned to consecutive time moments n=0,…,N.The map snapshotsMn are thus assigned to each of the temporary workspaces Pn. This is, however, difficult to be portrayed on the composite scheme given in Figure 1. Instead, there are only shaded planes illustrating this effect. A better portray of map snapshots is given in the following example (see Figure 2).Figure 2
Autonomous dynamics map of the process for the time momentt=0. Background color: the absolute value of the autonomous dynamics vector; small arrows: the direction of the autonomous dynamics vector; dotted lines: the orbits of the centers of attraction; solid lines: the borders of the forbidden zones; big arrows: the direction of the movement of the centers of attraction.Definition 2.11 (continuous-time dynamic system with discrete-time variability).
The temporary (potentially time-variant) system dynamics composed of the autonomous and forced dynamics can be defined in the computational moments as(2.9)x˙(tn)=x˙x(tn)+x˙u(tn)
which by the power of (2.3) and (2.4) allows for simple modeling of the system variability (which can be arbitrary induced by redefining the functions fx and fu).Note that the snapshotsMn of the map are included in (2.9) as x˙x, whereas the control part x˙u represents the forced component (2.4) of the complete state derivative x˙. Clearly, this means that not just the snapshot Mn forms the derivative of the state.Definition 2.12 (process computation).
The procedure of calculating the operational-point trajectory in the state-space for the assumed discrete-time variability model (which can also be treated as a kind of real-system approximation) consists of the following derivative-based prediction and numerical integration steps:(i)
takex˙(tn) and treat it as a local estimation of the system dynamics for the future time interval 〈tn,tn+1),(ii)
calculate the next state vector for the time momenttn+1 as
(2.10)x(tn+1)=x(tn)+hx˙(tn).Note that the rectangular integration scheme applied above for modeling the time-variability of dynamical systems is widely utilized in the procedures of computer simulation of variant or invariant continuous-time systems [7].Definition 2.13 (segmentation of the temporary operational subspaces).
Segmentation of the temporary operational subspacePn consists in dividing it, for the time moments tn, n=0,…,N, into a set 𝒮n of  S segments denoted as Φjn, for j=1,…,S, according to the following relation:
(2.11)Pn=⋃j=1S‍Φjn,
where each pair of segments for any discrete time n satisfies the following separability condition:
(2.12)(Φin-δΦin)∩(Φjn-δΦjn)=∅
for i,j=1,…,S;i≠j, where δ denotes the closure of the corresponding set. Clearly, the size of the segments is typically selected in accordance with some fidelity criteria of the applied discrete approximation of the original process.For the computational purposes we need a definition of a segment neighborhood and a generic concept of a distance (norm). Let us thus assume the common Euclidean norm∥·∥e and the following definition.Definition 2.14 (segment’s neighborhood).
Any segmentΦBn is in the neighborhood of another segment ΦAn, that is, ΦBn∈ℋ(ΦAn), if their geometrical centers xB and xA satisfy a simple neighboring condition:
(2.13)∥xB-xA∥e≤r,
where ∥·∥e means the Euclidean norm, and r is a fixed neighborhood size.Definition 2.15 (representatives).
Any variable used in optimization and assigned to a segment is said to be a representative.The autonomous dynamics map of Definition2.8, being the derivative of the vector function x(t), makes a principal source of candidates for a rational selection of the segment representatives. Note that in accordance to the common idea of domain and codomain of the transformation fx, from a practical implementation viewpoint, the map means a kind of doubling of the problem dimensioning with respect to Pn. This fact is merely represented in Figure 1 by the colors (or rather shades) applied to different temporary subspaces Pn. When the local evolution x˙ of some of the coordinates of the state vector is not relevant for the optimization considered, the co-domain can be reduced. At the same time the designer can introduce other representatives relevant for the optimization procedure and criterion.Definition 2.16 (predecessor and successor workspaces).
For any pair(Pn,Pn+1), Pn is called a predecessor and Pn+1 is referred to as a successor (on the time scale).Definition 2.17 (time-and-state operational space).
The time-and-state operational (sub)space, concerning all the time momentstn, n=0,…,N, is composed as follows:
(2.14)𝒫={P1,P2,…,Pn,…,PN}.
In spite of the above definition of the (global) time-and-state operational space, and according to the scheme depicted in Figure1, the principal optimization takes place in particular temporary workspaces Pn. Clearly, based on the procedure judgment shown in Figure 1, the time domain has to be treated separately.Definition 2.18 (state-space graph𝒢n).
A flow graph representing the arrangement of the segments of the temporary operational workspacePn, n=0,…,N, is said to be a state-space graph. According to the common idea [8] of flow graphs:
(2.15)𝒢n=(𝒲n,𝒟n,𝒦n),
we interpret its basic structure as a set 𝒲n of nodes, each of which is assigned a set of the values of the representatives of its corresponding segments. Only the edges connecting two nodes representing segments located in their mutual neighborhood form a set 𝒟n of the edges. Each edge d∈𝒟n is assigned a flow value k∈𝒦n.Definition 2.19 (time-and-state space toolgraph).
According to the definition of the time-and-state operational space, and based on the state-space graphs of all the workspacesPn, a time-and-state space toolgraph is composed by adding costless, directed edges connecting all the pairs of equivalent (in the subspace geometry) segments in the predecessor and successor workspaces.


### 2.3. Algorithm Description

Basically, the presented algorithm takes the above-defined principal notions (the time-and-state space𝒫 composed of the temporary operational workspaces Pn, the transition vector Λ, and the segmentation form) to minimize the adopted cost function J(E(u)). Optionally, we can determine the feasible control set 𝒰 and the forbidden zone Z. The latter can be used to represent some constraints portraying, for instance, faults previously detected in the system.Definition 2.20 (the toolgraph optimization algorithm).
Initialization with the input data: dynamic process description, operational time-and-state space, forbidden zones, initial and terminal trajectory points, and time horizon is as follows.(i)
Step  1. Segmentation.(ii)
Step  2. Determination of the representatives.(iii)
Step  3. Are the fidelity conditions met?
(1)
If yes, go to Step  4.(2)
If no, resize the segments, introduce new forbidden zones, and go to Step  1.(iv)
Step  4. Creation of the toolgraph.(v)
Step  5. T-graph search for the cheapest path connecting the initial and the terminal nodes.(vi)
Step  6. Has  the cheapest path been found?
(1)
If yes, go to Step  7.(2)
If no, terminate the algorithm.(vii)
Step  7. The reference trajectory synthesis (based on the path found in the previous step).As can be seen in Definition2.20, we start the procedure by segmenting the temporary subspaces Pn which results in a limited set of S segments, for n=0,…,N. Each segment of a set of representatives is assigned. Certainly, the size of the segments is chosen in accordance to an assumed discretization error, or fidelity index (discussed later on in the examples).Next, the state space graphs are formed. According to Definition2.18, each node from 𝒲n is associated with a certain segment of Pn. Its edges connect only nodes representing the segments neighboring in Pn. The edge flow values are determined according to the adopted function cost. Consequently, the time-and-state space toolgraph is built as a composition of the “consecutive” state-space graphs according to Definition 2.19.For practical optimization problems, any discrete optimization algorithm, which is suitable for finding the cheapest path between two selected nodes, can be utilized. If such a path exists, it represents a sequence of segments, which should be visited by the optimal trajectory. This sequence of segments implies a reference trajectory, which can be tracked by an executive controller in a procedure of a suboptimal fulfillment of the stated control task.


## 2.1. Problem Statement

Definition 2.1 (autonomous and forced dynamics).
A dynamical and time-variant process is described by the following function:(2.1)x˙(t)=f(x(t),u(t),t),(2.2)f:ℛl×ℛm×ℛ→ℛl,
where x=[x1,…,xl] is a state vector and u=[u1,…,um] is a control vector, and l,m∈𝒩. Let the autonomous dynamics be described by
(2.3)x˙x(t)=fx(x(t),t)
and let the forced dynamics be defined by a simple centroaffine transformation
(2.4)x˙u(t)=fu(x(t),u(t))=Fu(t),
where F is a fixed (time-invariant) l×m input control matrix, and the above contributors satisfy the following additivity condition:
(2.5)x˙(t)=x˙x(t)+x˙u(t),
where xx and xu are two, autonomous and forced state components.Definition 2.2 (operational subspaceP).
A subsetP=D1×D2×⋯×Dl, where Di, for i=1,…,l, represents a range along the axis xi, which is taken into account while seeking an optimal trajectory, is said to be an operational subspace.Definition 2.3 (forbidden zoneZ).
A subsetZ of P prohibited for operational points is referred to as a forbidden zone. This means that the sought optimal trajectory cannot enter it.Definition 2.4 (system trajectoryE(u)).
Any sequenceE(u) of consecutive states of a given dynamical system for a feasible control input u from a set 𝒰 is said to be a system’s trajectory or a trajectory of operational points in the system state space P∖Z.Definition 2.5 (transition vectorΛ).
A transition vectorΛ is an ordered set of two elements {x0,xk}(2.6)Λ={(x0,xk):x0,xk∈E(u)},
where x0 is the first element and xk is the last element of the sought (optimal) trajectory.Definition 2.6 (cost functionJ(E(u)) of the system and control trajectory).
Let𝒜 be a set of all trajectories E∈𝒜⊂P∖Z. Any real function of the class 𝒜→ℛ+ is said to be a cost function J(E(u)) of this trajectory.Definition 2.7 (optimal trajectoryE*).
LetΞ⊂𝒜 be the subset of all possible trajectories that perform the desired transition Λ (i.e., which start at x0 and terminate at xk). We say that a trajectory E* is (costly) optimal if it satisfies the following conditions:
(2.7)∀E∈ΞJ(E)≥J(E*),∀x∈E*x∈P∖Z.Definition 2.8 (autonomous dynamics mapM).
The imagex˙x of the whole space P in the transformation fx is said to be an autonomous dynamics map M(=x˙x).Our ultimate objective is to find the optimal and feasible controlu*(t) along with its corresponding state trajectory E*=E(u*(t))∈(P∖Z), which implements the complete transition Λ of the dynamical system (2.1) from its initial state x(t0)=x0 to the target state x(tk)=xk and at the same time minimizes the cost functional of the following form:
(2.8)J(E(u))=∫t0tk‍(∑i=1m‍βi|ui(t)|)dt,
where βi, i=1,…,m are nonnegative factors often used for weighting, and tk is a suitable transition-time interval resulting from the applied optimal control procedure.


## 2.2. Solution Method

For the dynamical systems described in Definition2.1, we consider its simple discrete-time variability of the computationally-convenient jump type allowing for the state-space system approximation scheme depicted in Figure 1.Figure 1
Continuous-time dynamic system with discrete-time variability.Leaving apart the initial conditions defined for a certaint0, let us consider consecutive continuous-time computational time moments tn, for n=1,…,N, with tn-tn-1=h=Δt, being a computation step (or a simulation insight), where N∈𝒩 is a discrete-time horizon under consideration.Definition 2.9 (temporary workspacePn).
At each of the above-mentioned computational moments, a separate temporary duplicatePn of the operational subspace P will be associated, as shown in Figure 1.Definition 2.10 (snapshotsMn of the dynamics map).
The multiplication of the subspaceP allows us to define separate snapshots Mn of the dynamics map M assigned to consecutive time moments n=0,…,N.The map snapshotsMn are thus assigned to each of the temporary workspaces Pn. This is, however, difficult to be portrayed on the composite scheme given in Figure 1. Instead, there are only shaded planes illustrating this effect. A better portray of map snapshots is given in the following example (see Figure 2).Figure 2
Autonomous dynamics map of the process for the time momentt=0. Background color: the absolute value of the autonomous dynamics vector; small arrows: the direction of the autonomous dynamics vector; dotted lines: the orbits of the centers of attraction; solid lines: the borders of the forbidden zones; big arrows: the direction of the movement of the centers of attraction.Definition 2.11 (continuous-time dynamic system with discrete-time variability).
The temporary (potentially time-variant) system dynamics composed of the autonomous and forced dynamics can be defined in the computational moments as(2.9)x˙(tn)=x˙x(tn)+x˙u(tn)
which by the power of (2.3) and (2.4) allows for simple modeling of the system variability (which can be arbitrary induced by redefining the functions fx and fu).Note that the snapshotsMn of the map are included in (2.9) as x˙x, whereas the control part x˙u represents the forced component (2.4) of the complete state derivative x˙. Clearly, this means that not just the snapshot Mn forms the derivative of the state.Definition 2.12 (process computation).
The procedure of calculating the operational-point trajectory in the state-space for the assumed discrete-time variability model (which can also be treated as a kind of real-system approximation) consists of the following derivative-based prediction and numerical integration steps:(i)
takex˙(tn) and treat it as a local estimation of the system dynamics for the future time interval 〈tn,tn+1),(ii)
calculate the next state vector for the time momenttn+1 as
(2.10)x(tn+1)=x(tn)+hx˙(tn).Note that the rectangular integration scheme applied above for modeling the time-variability of dynamical systems is widely utilized in the procedures of computer simulation of variant or invariant continuous-time systems [7].Definition 2.13 (segmentation of the temporary operational subspaces).
Segmentation of the temporary operational subspacePn consists in dividing it, for the time moments tn, n=0,…,N, into a set 𝒮n of  S segments denoted as Φjn, for j=1,…,S, according to the following relation:
(2.11)Pn=⋃j=1S‍Φjn,
where each pair of segments for any discrete time n satisfies the following separability condition:
(2.12)(Φin-δΦin)∩(Φjn-δΦjn)=∅
for i,j=1,…,S;i≠j, where δ denotes the closure of the corresponding set. Clearly, the size of the segments is typically selected in accordance with some fidelity criteria of the applied discrete approximation of the original process.For the computational purposes we need a definition of a segment neighborhood and a generic concept of a distance (norm). Let us thus assume the common Euclidean norm∥·∥e and the following definition.Definition 2.14 (segment’s neighborhood).
Any segmentΦBn is in the neighborhood of another segment ΦAn, that is, ΦBn∈ℋ(ΦAn), if their geometrical centers xB and xA satisfy a simple neighboring condition:
(2.13)∥xB-xA∥e≤r,
where ∥·∥e means the Euclidean norm, and r is a fixed neighborhood size.Definition 2.15 (representatives).
Any variable used in optimization and assigned to a segment is said to be a representative.The autonomous dynamics map of Definition2.8, being the derivative of the vector function x(t), makes a principal source of candidates for a rational selection of the segment representatives. Note that in accordance to the common idea of domain and codomain of the transformation fx, from a practical implementation viewpoint, the map means a kind of doubling of the problem dimensioning with respect to Pn. This fact is merely represented in Figure 1 by the colors (or rather shades) applied to different temporary subspaces Pn. When the local evolution x˙ of some of the coordinates of the state vector is not relevant for the optimization considered, the co-domain can be reduced. At the same time the designer can introduce other representatives relevant for the optimization procedure and criterion.Definition 2.16 (predecessor and successor workspaces).
For any pair(Pn,Pn+1), Pn is called a predecessor and Pn+1 is referred to as a successor (on the time scale).Definition 2.17 (time-and-state operational space).
The time-and-state operational (sub)space, concerning all the time momentstn, n=0,…,N, is composed as follows:
(2.14)𝒫={P1,P2,…,Pn,…,PN}.
In spite of the above definition of the (global) time-and-state operational space, and according to the scheme depicted in Figure1, the principal optimization takes place in particular temporary workspaces Pn. Clearly, based on the procedure judgment shown in Figure 1, the time domain has to be treated separately.Definition 2.18 (state-space graph𝒢n).
A flow graph representing the arrangement of the segments of the temporary operational workspacePn, n=0,…,N, is said to be a state-space graph. According to the common idea [8] of flow graphs:
(2.15)𝒢n=(𝒲n,𝒟n,𝒦n),
we interpret its basic structure as a set 𝒲n of nodes, each of which is assigned a set of the values of the representatives of its corresponding segments. Only the edges connecting two nodes representing segments located in their mutual neighborhood form a set 𝒟n of the edges. Each edge d∈𝒟n is assigned a flow value k∈𝒦n.Definition 2.19 (time-and-state space toolgraph).
According to the definition of the time-and-state operational space, and based on the state-space graphs of all the workspacesPn, a time-and-state space toolgraph is composed by adding costless, directed edges connecting all the pairs of equivalent (in the subspace geometry) segments in the predecessor and successor workspaces.


## 2.3. Algorithm Description

Basically, the presented algorithm takes the above-defined principal notions (the time-and-state space𝒫 composed of the temporary operational workspaces Pn, the transition vector Λ, and the segmentation form) to minimize the adopted cost function J(E(u)). Optionally, we can determine the feasible control set 𝒰 and the forbidden zone Z. The latter can be used to represent some constraints portraying, for instance, faults previously detected in the system.Definition 2.20 (the toolgraph optimization algorithm).
Initialization with the input data: dynamic process description, operational time-and-state space, forbidden zones, initial and terminal trajectory points, and time horizon is as follows.(i)
Step  1. Segmentation.(ii)
Step  2. Determination of the representatives.(iii)
Step  3. Are the fidelity conditions met?
(1)
If yes, go to Step  4.(2)
If no, resize the segments, introduce new forbidden zones, and go to Step  1.(iv)
Step  4. Creation of the toolgraph.(v)
Step  5. T-graph search for the cheapest path connecting the initial and the terminal nodes.(vi)
Step  6. Has  the cheapest path been found?
(1)
If yes, go to Step  7.(2)
If no, terminate the algorithm.(vii)
Step  7. The reference trajectory synthesis (based on the path found in the previous step).As can be seen in Definition2.20, we start the procedure by segmenting the temporary subspaces Pn which results in a limited set of S segments, for n=0,…,N. Each segment of a set of representatives is assigned. Certainly, the size of the segments is chosen in accordance to an assumed discretization error, or fidelity index (discussed later on in the examples).Next, the state space graphs are formed. According to Definition2.18, each node from 𝒲n is associated with a certain segment of Pn. Its edges connect only nodes representing the segments neighboring in Pn. The edge flow values are determined according to the adopted function cost. Consequently, the time-and-state space toolgraph is built as a composition of the “consecutive” state-space graphs according to Definition 2.19.For practical optimization problems, any discrete optimization algorithm, which is suitable for finding the cheapest path between two selected nodes, can be utilized. If such a path exists, it represents a sequence of segments, which should be visited by the optimal trajectory. This sequence of segments implies a reference trajectory, which can be tracked by an executive controller in a procedure of a suboptimal fulfillment of the stated control task.


## 3. Optimization Examples

This section provides an illustrative example which utilizes the above-mentioned formal definitions and the optimization algorithm. The minimized criterion is given by (2.8) with the restriction that the resulting trajectory is piecewise linear in the considered time-and-state space. The example concerns a two-dimensional time-variant dynamical process described in the following subsection.


### 3.1. Process Description

The exemplary process is described by the following set of equations:(3.1)x˙1(t)=GMAx1dA(x1(t),x2(t))2-GMB(x1(t)-x1B(t))dB(x1(t),x2(t))2-GMC(x1(t)-x1C(t))dC(x1(t),x2(t))2-GMD(x1(t)-x1D(t))dD(x1(t),x2(t))2,x˙2(t)=GMAx2dA(x1(t),x2(t))2-GMB(x2(t)-x2B(t))dB(x1(t),x2(t))2-GMC(x2(t)-x2C(t))dC(x1(t),x2(t))2-GMD(x2(t)-x2D(t))dD(x1(t),x2(t))2,
where G=6.673×10-3, MA=7.9891×102, MB=4.8685×102, MC=8.9736×102, and MD=7.3476×102 are the process constant parameters. The auxiliary model variables and functions given above are defined as follows:
(3.2)dA(x1(t),x2(t))=(x1(t))2+(x2(t))2,dB(x1(t),x2(t))=(x1(t)-x1B(t))2+(x2(t)-x2B(t))2,dC(x1(t),x2(t))=(x1(t)-x1C(t))2+(x2(t)-x2C(t))2,dD(x1(t),x2(t))=(x1(t)-x1D(t))2+(x2(t)-x2D(t))2,
where
(3.3)x1B(t)=RBsinαB(t),x2B(t)=RBcosαB(t),x1C(t)=RCsinαC(t),x2C(t)=RCcosαC(t),x1D(t)=RDsinαD(t)+RCsinαC(t),x2D(t)=RDcosαD(t)+RCcosαC(t),RB=50,RC=120,RD=40
and αB, αC, αD are the internal process state variables characterized by derivatives
(3.4)α˙B(t)=ΩB,α˙C(t)=ΩC,α˙D(t)=ΩD
of the following constant values: ΩB=65.306×10-3, ΩC=17.56×10-3,  ΩD=96.728×10-3.The autonomous dynamics map of the described process is shown in Figure2, where certain three (B,C,D) centers of attraction can be seen, which circulating with the angular velocities (ΩB, ΩC, ΩD, resp.) around the center (A) are the source of process time-variability. Additionally, the forbidden zones ZA, ZB, ZC, ZD are defined in the forms of circular areas, which have the radii rA=10, rB=20, rC=15, and rD=10, respectively. The midpoints of the forbidden zones cover the midpoints of the respective centers of “attraction.”


### 3.2. Algorithm Implementation

In theory, an infinite time horizon of optimization is a simple and trivial solution for the time horizon setting. In computer applications, we prefer the concept of “sufficiently” large horizons, which is determined according to our prior knowledge on the optimization process and the limitations of available computing power, whereas the time step size was chosen empirically so as to obtain the desired “granularity” of the solution.According to Definition2.13, the segmentation procedure consists of dividing the operational subspaces Pn, n=0,…,N, being, in this case, hypercubes described by D1×D2=(-250,250)×(-250,250), into the set of segments of diameters (rx1=10,rx2=10), all located along the time axis at a step size Δt=10. For each segment the representative values are calculated as the averages of the autonomous dynamics computed for particular directions:
(3.5)FAVRΦj,in=1Np∑k=1Np‍FDYNi,kn,
where FDYNi,kn denotes the kth sample of the coordinate xi of the autonomous dynamics map for the discrete time moments n, and Np=25 denotes the number of samples per segment in this following examples.The fidelity of the applied discrete representation for the segmentΦj,in is assessed as
(3.6)eΦj,in=∑k=1Np‍|FDYNi,kn-FAVRΦj,in|.IfeΦj,in exceeds some given threshold, which can be an empirically estimated value, we assume that the segment Φj,in along with its representatives constitutes an infeasible representation of the autonomous dynamics in the area of the time-and-state space occupied by this segment. In such cases either the size of the segment is changed or the segment is included to the forbidden zone (as an area of high dynamics).Having defined the neighborhood of the segments, the construction of the state-space graphs and time-and-state space toolgraph comes according to Definitions2.18 and 2.19, respectively. The flow values are determined as the cost of driving the operational point between each two neighboring segments according to (2.8).The last stage consists in searching the above-defined time-and-state space toolgraph for the cheapest path connecting the initial and the terminal node. We can use, for instance, the Dijkstra’s algorithm [9] for this purpose.


### 3.3. Problem 1

The first problem concerns finding a trajectory of the operational point, which minimizes (2.8), with the initial location at (x1o=-250,x2o=-250,t0=0) and the terminal point (x1k,x2k,tk) satisfying the following condition:
(3.7)(x1k-x1D)2+(x2k-x2D)2≤1.5rd,(x1k,x2k)∉ZD,tk≤tN,
where the assumed time horizon is tN=2000. The results are presented in Figure 3. The total cost indicator is equal to J1=43.454×10-2. It can be clearly seen that the optimal strategy consists in keeping the operational point at the initial location till a certain time moment (t*=1405) and then to start the transition process.Figure 3
Resulting optimal trajectory for Problem 1 (solid black line); the forbidden zones of the respective centers of attractions:A, B, C, D (colored areas).


### 3.4. Problem 2

In this case the time horizon was shortened totN=1600, which resulted in a higher value J2=52.494×10-2 of the total cost indicator as compared to the one of Problem 1. The remaining parameters are the same as in Problem 1. The results are depicted in Figure 4.Figure 4
Optimal trajectory for Problem 2 (solid black line); the forbidden zones of the respective centers of attractions:A, B, C, D (colored areas).


### 3.5. Problem 3

Further shortening of the time horizontN=1000 leads to a greater control cost J3=150.744×10-2 with the trajectory of the operational point depicted in Figure 5.Figure 5
Optimal trajectory for Problem 3 (solid black line); the forbidden zones of the respective centers of attractions:A, B, C, D (colored areas).


### 3.6. Problem 4

This time the problem of driving the operational point from the initial condition(x1o=0,x2o=0,t0=0) to the terminal point (x1k>180,x2k>180,tk≤tN) is considered. The time horizon is tN=1000. The resulting total cost is equal to J4=275.368×10-2. The designed trajectory is presented in Figure 6.Figure 6
Optimal trajectory for Problem 4 (solid black line); the forbidden zones of the respective centers of attractions:A, B, C, D (colored areas).


### 3.7. Problem 5

This case is a variation of Problem 4 with the time horizon lengthened totN=4000, which resulted in a reduction of the total cost J5=202.905×10-2. The computed trajectory is presented in Figure 7.Figure 7
Optimal trajectory for Problem 5 (solid black line); the forbidden zones of the respective centers of attractions:A, B, C, D (grayed areas).


### 3.8. Problem 6

The last example presents the design of the optimal trajectory with the initial point at(x1o=2.5,x2o=-17.5,t0=0) and the terminal point satisfying (3.7) with the time horizon tN=4000. The total cost J5=24.752×10-2. The gained results are depicted in Figure 8.Figure 8
Optimal trajectory for Problem 6 (solid black line); the forbidden zones of the respective centers of attractions:A, B, C, D (colored areas).


## 3.1. Process Description

The exemplary process is described by the following set of equations:(3.1)x˙1(t)=GMAx1dA(x1(t),x2(t))2-GMB(x1(t)-x1B(t))dB(x1(t),x2(t))2-GMC(x1(t)-x1C(t))dC(x1(t),x2(t))2-GMD(x1(t)-x1D(t))dD(x1(t),x2(t))2,x˙2(t)=GMAx2dA(x1(t),x2(t))2-GMB(x2(t)-x2B(t))dB(x1(t),x2(t))2-GMC(x2(t)-x2C(t))dC(x1(t),x2(t))2-GMD(x2(t)-x2D(t))dD(x1(t),x2(t))2,
where G=6.673×10-3, MA=7.9891×102, MB=4.8685×102, MC=8.9736×102, and MD=7.3476×102 are the process constant parameters. The auxiliary model variables and functions given above are defined as follows:
(3.2)dA(x1(t),x2(t))=(x1(t))2+(x2(t))2,dB(x1(t),x2(t))=(x1(t)-x1B(t))2+(x2(t)-x2B(t))2,dC(x1(t),x2(t))=(x1(t)-x1C(t))2+(x2(t)-x2C(t))2,dD(x1(t),x2(t))=(x1(t)-x1D(t))2+(x2(t)-x2D(t))2,
where
(3.3)x1B(t)=RBsinαB(t),x2B(t)=RBcosαB(t),x1C(t)=RCsinαC(t),x2C(t)=RCcosαC(t),x1D(t)=RDsinαD(t)+RCsinαC(t),x2D(t)=RDcosαD(t)+RCcosαC(t),RB=50,RC=120,RD=40
and αB, αC, αD are the internal process state variables characterized by derivatives
(3.4)α˙B(t)=ΩB,α˙C(t)=ΩC,α˙D(t)=ΩD
of the following constant values: ΩB=65.306×10-3, ΩC=17.56×10-3,  ΩD=96.728×10-3.The autonomous dynamics map of the described process is shown in Figure2, where certain three (B,C,D) centers of attraction can be seen, which circulating with the angular velocities (ΩB, ΩC, ΩD, resp.) around the center (A) are the source of process time-variability. Additionally, the forbidden zones ZA, ZB, ZC, ZD are defined in the forms of circular areas, which have the radii rA=10, rB=20, rC=15, and rD=10, respectively. The midpoints of the forbidden zones cover the midpoints of the respective centers of “attraction.”


## 3.2. Algorithm Implementation

In theory, an infinite time horizon of optimization is a simple and trivial solution for the time horizon setting. In computer applications, we prefer the concept of “sufficiently” large horizons, which is determined according to our prior knowledge on the optimization process and the limitations of available computing power, whereas the time step size was chosen empirically so as to obtain the desired “granularity” of the solution.According to Definition2.13, the segmentation procedure consists of dividing the operational subspaces Pn, n=0,…,N, being, in this case, hypercubes described by D1×D2=(-250,250)×(-250,250), into the set of segments of diameters (rx1=10,rx2=10), all located along the time axis at a step size Δt=10. For each segment the representative values are calculated as the averages of the autonomous dynamics computed for particular directions:
(3.5)FAVRΦj,in=1Np∑k=1Np‍FDYNi,kn,
where FDYNi,kn denotes the kth sample of the coordinate xi of the autonomous dynamics map for the discrete time moments n, and Np=25 denotes the number of samples per segment in this following examples.The fidelity of the applied discrete representation for the segmentΦj,in is assessed as
(3.6)eΦj,in=∑k=1Np‍|FDYNi,kn-FAVRΦj,in|.IfeΦj,in exceeds some given threshold, which can be an empirically estimated value, we assume that the segment Φj,in along with its representatives constitutes an infeasible representation of the autonomous dynamics in the area of the time-and-state space occupied by this segment. In such cases either the size of the segment is changed or the segment is included to the forbidden zone (as an area of high dynamics).Having defined the neighborhood of the segments, the construction of the state-space graphs and time-and-state space toolgraph comes according to Definitions2.18 and 2.19, respectively. The flow values are determined as the cost of driving the operational point between each two neighboring segments according to (2.8).The last stage consists in searching the above-defined time-and-state space toolgraph for the cheapest path connecting the initial and the terminal node. We can use, for instance, the Dijkstra’s algorithm [9] for this purpose.


## 3.3. Problem 1

The first problem concerns finding a trajectory of the operational point, which minimizes (2.8), with the initial location at (x1o=-250,x2o=-250,t0=0) and the terminal point (x1k,x2k,tk) satisfying the following condition:
(3.7)(x1k-x1D)2+(x2k-x2D)2≤1.5rd,(x1k,x2k)∉ZD,tk≤tN,
where the assumed time horizon is tN=2000. The results are presented in Figure 3. The total cost indicator is equal to J1=43.454×10-2. It can be clearly seen that the optimal strategy consists in keeping the operational point at the initial location till a certain time moment (t*=1405) and then to start the transition process.Figure 3
Resulting optimal trajectory for Problem 1 (solid black line); the forbidden zones of the respective centers of attractions:A, B, C, D (colored areas).


## 3.4. Problem 2

In this case the time horizon was shortened totN=1600, which resulted in a higher value J2=52.494×10-2 of the total cost indicator as compared to the one of Problem 1. The remaining parameters are the same as in Problem 1. The results are depicted in Figure 4.Figure 4
Optimal trajectory for Problem 2 (solid black line); the forbidden zones of the respective centers of attractions:A, B, C, D (colored areas).


## 3.5. Problem 3

Further shortening of the time horizontN=1000 leads to a greater control cost J3=150.744×10-2 with the trajectory of the operational point depicted in Figure 5.Figure 5
Optimal trajectory for Problem 3 (solid black line); the forbidden zones of the respective centers of attractions:A, B, C, D (colored areas).


## 3.6. Problem 4

This time the problem of driving the operational point from the initial condition(x1o=0,x2o=0,t0=0) to the terminal point (x1k>180,x2k>180,tk≤tN) is considered. The time horizon is tN=1000. The resulting total cost is equal to J4=275.368×10-2. The designed trajectory is presented in Figure 6.Figure 6
Optimal trajectory for Problem 4 (solid black line); the forbidden zones of the respective centers of attractions:A, B, C, D (colored areas).


## 3.7. Problem 5

This case is a variation of Problem 4 with the time horizon lengthened totN=4000, which resulted in a reduction of the total cost J5=202.905×10-2. The computed trajectory is presented in Figure 7.Figure 7
Optimal trajectory for Problem 5 (solid black line); the forbidden zones of the respective centers of attractions:A, B, C, D (grayed areas).


## 3.8. Problem 6

The last example presents the design of the optimal trajectory with the initial point at(x1o=2.5,x2o=-17.5,t0=0) and the terminal point satisfying (3.7) with the time horizon tN=4000. The total cost J5=24.752×10-2. The gained results are depicted in Figure 8.Figure 8
Optimal trajectory for Problem 6 (solid black line); the forbidden zones of the respective centers of attractions:A, B, C, D (colored areas).


## 4. Conclusions

This work presents an effective concept of designing optimal and feasible control strategies for time-varying dynamical processes based on discrete optimization approach and its solutions.After segmenting the composed time-and-state workspace of the process dynamics, a time-and-state space toolgraph (made of temporary state-space graphs) is constructed that represents all the pertinent properties of the autonomous dynamics of the considered dynamical process. Additional assumptions or restrictions concerning arbitrary forbidden zones for the operational points can be easily taken into account. Numerical computations are performed by means of any discrete optimization algorithm which searches for the cheapest path connecting the node initial and terminal nodes. The presented work has been implemented based on the basic form of the Dijkstra’s algorithm, whose performance can be improved by applying various modifications given in the literature, for instance, in [10].The resulting cheapest path describes an optimal piecewise-linear trajectory of the operational points in the considered time-and-state space. Note that such a solution may need some kind of smoothing for most real-time process cases.A clear suboptimal property of the proposed approach results from the  fact that the computed reference trajectory can only pass through the centers of the segments. In this context the size of segments has a key influence on the quality of the ultimate optimal solution being sought.Thus, in general, this interdisciplinary paper concerns research into formal optimization methods and tools. Namely, it shows how to constructively predefine the control optimisation problem or to convert complex optimal control problems for nonlinear time-variant objects with constraints into a graph structure ready to be processed by any of the discrete optimization methods, in particular, by the standard numerical methods of graph optimisation. In the light of theory application to problems arising in engineering, this method can be used for on- or off-line numerically solved practical optimal control or path-finding problems for time-varying objects or environments.Though this paper is an investigative proposition of a new numerical optimization method, an abstract but nontrivial and imaginative example of several variations has also been given to illustrate the effectiveness of the proposed methodology. The example presented inclusive of a system of variant attraction centers is a rough analog to the solar system. Other business and industry applications are to be developed in order to give a wider proof of the advantages of the proposed method in a confrontation with the scientific challenges of the real world. We also hope that this paper should encourage scientists and engineers to apply both this approach and the discrete optimisation methods.Complexity of the presented method principally results from the complexity of the algorithm applied for searching the toolgraph. Some clues concerning the possibilities of improving the graph search process can be found in [11, 12], for instance. It is clear that the difficulty of the proposed method is shaped mostly by the numerical complexity of searching the T-graph. In a straightforward tactic one makes use of the Dijkstra algorithm (O(n2)). Other, mostly heuristic enhancements (like the renowned algorithm A*) can also be employed in order to improve the efficiency of the proposed numerical optimisation procedure.In spite of the presented feasibility of the approach, a general problem of the so-called curse of dimensionality has always to be taken into account in all research of this type. This means that there are limits which render such approaches infeasible in some problems of practical interest. Nevertheless, nowadays, in view of the persistent progress in computer technology, finding constructive methods of solving problems appears to be most substantial in science, technology, and engineering.



---

*Source: 142060-2012-10-31.xml*