# Approximate Solution of a Reduced-Type Index-k Hessenberg Differential-Algebraic Control System

Authors: Ghazwa F. Abd; Radhi A. Zaboon
Journal: Journal of Applied Mathematics
 (2021)

---

## Abstract

This study focuses on developing an efficient and easily implemented novel technique to solve the index-k Hessenberg differential-algebraic equation (DAE) with input control. The implicit function theorem is first applied to solve the algebraic constraints of having unknown state differential variables to form a reduced state-space representation of an ordinary differential (control) system defined on smooth manifold with consistent initial conditions. The variational formulation is then developed for the reduced problem. A solution of the reduced problem is proven to be the critical point of the variational formulation, and the critical points of the variational formulation are the solutions of the reduced problem on the manifold. The approximate analytical solution of the equivalent variational formulation is represented as a finite number of basis functions with unknown parameters on a suitable separable Hilbert setting solution space. The unknown coefficients of the solution are obtained by solving a linear algebraic system. The different index problems of linear Hessenberg differential-algebraic control systems are approximately solved using this approach with comparisons. The numerical results reveal the good efficiency and accuracy of the proposed method. This technique is applicable for a large number of applications like linear quadratic optima, control problems, and constrained mechanical systems.


---

## Body

## 1. Introduction

The differential-algebraic equation system is often used to model many real-life application problems arising in electrical networks, optimal control systems, and constrained multibody mechanical systems [1–5]. The index property is often used to characterize a differential-algebraic system, and it has been given several definitions [3, 4, 6]. The differentiation index is the most commonly used and represents the minimum number of times that the algebraic constraint of DAEs must be differentiated with respect to time to obtain an equivalent ordinary differential equation [4]. High-index DAEs are more challenging to solve than low-index DAEs [7]. Therefore, index reduction transformation techniques have been applied to transform higher-index DAEs into an index of zero or one before solving them [8]. Some numerical methods have been developed for solving special structure, low-index DAEs [9–14].Hessenberg DAEs are a very important class of general DAEs with a special structure and a wide variety of scientific and engineering applications [15, 16]. The homotopy analysis method has been effectively applied for some classes to approximate the analytical solution of DAE classes [8]. A numerical method for solving index 2 and index 3 Hessenberg DAEs is presented in [17–21]. Most of these methods are based on the index reduction method, which transforms a high DAE index into a low DAE index. A novel technique for solving nonlinear high-index DAEs using the Adomian decomposition method is proposed in [22]. Reliable techniques for solving high-index Hessenberg DAEs are required.Most of the previous methods were applied to the lower-order index problem (index less than 3) [23, 24]. The motivation of this work is to develop an easily implemented and efficient theoretical method to solve the higher-index Hessenberg control problem. The approximate analytical solution is easily obtained using the following steps: First, the implicit function theorem is applied to the algebraic constraints of the given DAEs to obtain a reduced state-space differential equation defined on a smooth manifold. Then, a variational formulation equivalent to the reduced problem is developed so that both problems have the same solution.Finally, the analytical solution of the variational formulation is approximated on a separable complete sitting space using a countable basis function representation with unknown variables, which are obtained as a solution of the solvable linear algebraic system. By solving this system, the approximate analytical solution can be easily obtained. This approach is efficiently applied to the high-index Hessenberg problem.


## 2. Problem Formulation of Index-k Hessenberg DAEs

One can consider the following index-k linear Hessenberg DAEs with control input [10]:
(1)I0..00I.....I.....I.0...0ẏ1⋮ẏk−1ẏk+A11..A1,k−1A1,kA21..A2,k−100.........0.0Ak,k−10y1⋮yk−1yk=B1⋮Bk−1Bku+f1⋮fk−1fk,where y1,⋯,ykT denotes independent state vectors; Aij;i,j=1,..,kare the constant matrices of appropriate dimensions, and Aij does not need to be a square matrix or invertible; Ak,k−1Ak−1,k−2⋯A1,k≜c is a square invertible matrix, where the index is the differentiation index (see [10, 16, 18]); Bj,j=1,..,k are appropriate constant vectors, and u is the given single control input with u∈ck−1T;R≜∆u the class of admissible control; T=t0,tf,t0<tf; fj;j=1,..,k denotes differentiable functions; fj∈ck−1T;R.The DAEs with the givenu∈∆u are linear time-invariant Hessenberg DAEs of size k, which are solvable for u∈∆u and have an index of k. (1) can be rewritten as follows:
(2)ẏ1+A11y1+⋯+A1,kyk=B1u+f1,ẏ2+A21y1+⋯+A2,k−1yk−1=B2u+f2,⋮ẏk−2+Ak−2,k−3yk−3+Ak−2,k−2yk−2+Ak−2,k−1yk−1=Bk−2u+fk−2,ẏk−1+Ak−1,k−2yk−2+Ak−1,k−1yk−1=Bk−1u+fk−1,Ak,k−1yk−1=Bku+fk.The square invertible matrixAk,k−1Ak−1,k−2⋯A1,k represents the Jacobean condition; thus, the implicit function theorem [25–27] can be applied to estimate the unknown state yk from the last algebraic equation using the differential index property. The algebraic constraint can be differentiated using Equation (2) with respect to t:
(3)Ak,k−1ẏk−1=Bku̇+ḟk.Substituting the differential equationẏk−1 from Equation (2) into Equation (3) results in the following:
(4)−Bku̇−ḟk+Ak,k−1Bk−1u+fk−1−Ak−1,k−2yk−2−Ak−1,k−1yk−1=0.Equation (4) now performs a new constraint conjugated with the given differential-algebraic system in Equation (2). Moreover, the state variable yk does not appear in Equation (4). Therefore, we can differentiate Equation (4) with respect to t as follows:
(5)−Bku¨−f¨k+Ak,k−1Bk−1u̇+ḟk−1−Ak−1,k−2ẏk−2−Ak−1,k−1ẏk−1=0.Next, one can substitute the differential equationsẏk−2 and ẏk−1 obtained from Equation (2):
(6)−Bku¨−f¨k+Ak,k−1Bk−1u̇+Ak,k−1ḟk−1−Ak,k−1Ak−1,k−2ẏk−2−Ak,k−1Ak−1,k−1ẏk−1=0⟹−Bku¨+Ak,k−1Bk−1u̇−Ak,k−1Ak−1,k−2Bk−2+Ak,k−1Ak−1,k−1Bk−1u−f¨k+Ak,k−1ḟk−1−Ak,k−1Ak−1,k−2fk−2−Ak,k−1Ak−1,k−1fk−1+Ak,k−1Ak−1,k−2Ak−2,k−3yk−3+Ak,k−1Ak−1,k−2Ak−2,k−2+Ak,k−1Ak−1,k−1Ak−1,k−2yk−2+Ak,k−1Ak−1,k−2Ak−2,k−1+Ak,k−1Ak−1,k−1Ak−1,k−1yk−1=0.Equation (6) represents another algebraic constraint added to system (2). The derivation is continued with respect to t until the index-k condition (Ak,k−1Ak−1,k−2⋯A1k) is satisfied to estimate the state variable yk:
(7)−Bkuk−1+Ak,k−1Bk−1uk−2−⋯−1kAk,k−1Ak−1,k−1⋯A21B1u−fkk−1+Ak,k−1fk−1k−2−⋯−1kAk,k−1Ak−1,k−1⋯A21f1−Ak,k−1Ak−1,k−2⋯A1kyk−s1yk−1−s2yk−2−⋯−sk−1y1=0.The matrixAk,k−1Ak−1,k−2⋯A1k≜c is nonsingular; thus,
(8)s1=Ak,k−1Ak−1,k−2Ak−2,k−3Ak−3,k−1+Ak,k−1Ak−1,k−2Ak−2,k−2+Ak,k−1Ak−1,k−1Ak−1,k−2Ak−2,k−1+Ak,k−1Ak−1,k−2Ak−2,k−1+Ak,k−1Ak−1,k−1Ak−1,k−1Ak−1,k−1,s2=Ak,k−1Ak−1,k−2Ak−2,k−3Ak−3,k−2+Ak,k−1Ak−1,k−2Ak−2,k−2+Ak,k−1Ak−1,k−1Ak−1,k−2Ak−2,k−2+Ak,k−1Ak−1,k−2Ak−2,k−1+Ak,k−1Ak−1,k−1Ak−1,k−1Ak−1,k−2⋯.Then,(9)yk=c−1−1k−1Bkuk−1+−1kAk,k−1Bk−1uk−2−⋯−1k−1Ak,k−1Ak−1,k−1⋯A21B1u+−1k−1fkk−1+−1kAk,k−1fk−1k−2−⋯−1k−1Ak,k−1Ak−1,k−1⋯A21f1−⋯−s1yk−1−s2yk−2−⋯−sk−1y1=Ly1,y2,⋯,yk−1+Lu,u̇,⋯,uk−1,fj,ḟj,⋯,fjk−1,j=1,..,k,where Ly1,y2,⋯,yk−1=−s1yk−1−s2yk−2−⋯−sk−1y1,
(10)Lu,u̇,⋯,uk−1,fj,ḟj,⋯,fjk−1=−1k−1Bkuk−1+−1kAk,k−1Bk−1uk−2−⋯−1k−1Ak,k−1Ak−1,k−1⋯A21B1u+−1k−1fkk−1+−1kAk,k−1fk−1k−2−⋯−1k−1Ak,k−1Ak−1,k−1⋯A21f1−⋯.System (1) is an index-k differential-algebraic system of the Hessenberg type. The initial condition should be chosen consistently, such that all the algebraic constraints in Equations (3), (4), and (6) are satisfied. The cost is calculated using the differentiation index to obtain an explicit linear ordinary differential equation defined on a manifold with consistent initial conditions that the time-dependent ut,fj,j=1,..,k is also differentiated k−1 times, influencing the solution with u,u̇,⋯,uk−1,fj,ḟj,⋯,fjk−1 as a system of input perturbation.From the additionalk−1 hidden constraint Equations (4) and (6) and yk=Ly1,y2,⋯,yk−1+Lu,u̇,⋯,uk−1,fj,ḟj,⋯,fjk−1, the reduced-type index-k Hessenberg DAE (2) with control u,u̇,⋯,uk−1 over the class of ck−1T;Rtransforms into the following:
(11)ẏ1+A11y1+⋯+A1,kLy1,y2,⋯,yk−1=B1u+f1+Lu~,f~,ẏ2+A21y1+⋯+A2,k−1yk−1=B2u+f2,⋮Ak,k−1yk−1=Bku+fk,Ak−1,k−2yk−2+Ak−1,k−1yk−1=−Bku̇+Ak,k−1Bk−1u−ḟk+Ak,k−1fk−1,⋮s2yk−1+s2yk−2+⋯+sk−2y1=Bk−1uk−2−⋯−1k−1Ak,k−1Ak−1,k−1⋯A21B1+fk−1k−2−⋯−1k−1Ak,k−1Ak−1,k−1⋯A21f1,Ak,k−1yk−1=Bku+fkt=t0,Ak−1,k−2yk−2+Ak−1,k−1yk−1=−Bku̇+Ak,k−1Bk−1u−ḟk+Ak,k−1fk−1t=t0,⋮s2yk−1+s2yk−2+⋯+sk−2y1=Bk−1uk−2−⋯−1k−1Ak,k−1Ak−1,k−1⋯A21B1u+fk−1k−2−⋯−1k−1Ak,k−1Ak−1,k−1⋯A21f1t=t0,where Lu~,f~=Lu,u̇,⋯,uk−1,fj,ḟj,⋯,fjk−1.


## 3. Solving the Reduced-Type Index-k Hessenberg DAE Using Variational Formulation

To solve the problem in Equation (11) using the variational formulation, one can define the operator L:DL⊂H⟶RL⊂H; with DL¯=H, where H is a suitable Hilbert space (e.g., H≜L2T;T=t0,tf,t0<tf) by the following:
(12)Ly=Ft;u~,f~,where the operator L is a linear differential and algebraic operator defined as follows:
(13)Ly=I0.00I....I.0.000.00⋮dy1dt⋮dyk−1dtdykdt+A11..A1,k−1A1A21..A2,k−100......0..0.0Ak,k−100.A~2A~30⋮y1⋮yk−1yk,(14)y=y1,⋯,yk,A~1=Ak,k−1Ly1,y2,⋯,yk−1,A~2=Ak,k−1Ak−1,k−1,A~3=−Ak,k−1Ak−1,k−2,(15)Ft;u~,f~≜Ft;u,u̇,u¨,⋯uk−1,fj,ḟj,⋯,fjk−1.For allu∈ck−1T;R≜∆u, the class of admissible control, one can define
(16)ω0=yt0=y1t0,..,ykt0T∈RnAk,k−1yk−1t0=Bkut0andAk−1,k−2yk−2t0+Ak−1,k−1yk−1t0=−Bku̇t0+Ak,k−1Bk−1ut0−ḟkt0+Ak,k−1fk−1t0⋮s2yk−1t0+s2yk−2t0+⋯+sk−2y1t0=Bk−1uk−2t0−⋯−1k−1Ak,k−1Ak−1,k−1⋯A21B1ut0+fk−1k−2t0−⋯−1k−1Ak,k−1Ak−1,k−1⋯A21f1t0andDL=yt=y1t,..,yktT∈Rk∣y∈L2T;Rk,yt0∈ω0.This system of linearly invariant reduced-type index Hessenberg DAEs satisfies the DAEs almost everywhere inH. The solution is assumed to be an absolutely continuous function satisfying (11) almost everywhere over the setting class C1T;Rk. In addition, yk is assumed to be continuously differentiable, satisfying Equation (11).


## 4. Variational Formulation Algorithm

The following algorithm is proposed to simplify the discussion of the presented approach:Step 1.
One can consider the index-k Hessenberg control DAEs given in Equation (1).Step 2.
The equivalent reduced-type index-k Hessenberg DAEs can be derived from Equation (11).Step 3.
Because of the presentd/dt operator, operator L defined in Equation (13) is not symmetric with respect to the usual bilinear form y,y~=∫yTt.y~tdt,y=y1t,..,ykt,y~=y~1t,⋯,y~kt. This bilinear form should be chosen to satisfy the nondegenerate condition on the domain and the range of the linear operator L (i.e., for every ∈DL,y,y~=0, y~=0, and for every y~∈DL,y,y~=0,y=0. Thus, no usual variational formulation exists unless one can redefine either the linear operator or the usual bilinear (inner product) form using the form y,y~≜y,Ly~,y,y~∈DL, where (DL⊂HL⊂H with yL=y,Ly, HL,dy,y~ is the complete Hilbert space, andHL is assumed to be a separable Hilbert space. Then, L is a symmetric linear operator (see [24–26]). If L is positive (positive definite) on DL¯=H, then minJy is defined; otherwise, the critical points of Jω0y should be obtained, where the augmented functional Jω0y over the class of the consistent initial conditions and u∈ck−1T;R is as follows:
(17)Jω0y=12L1T,L1−F1T,L1+⋯+12LkT,Lk−FkT,Lk+12Lk+1T,Lk+1−Fk+1T,Lk+1+12Lk+2T,Lk+2−Fk+2T,Lk+2+⋯+12Lk+k−1T,Lk+k−1−Fk+k−1T,Lk+k−1+12LkTt0,Lkt0−FkTt0,Lkt0+12Lk+1Tt0,Lk+1t0−Fk+1Tt0,Lk+1t0+⋯+12Lk+k−1Tt0,Lk+k−1t0−Fk+k−1Tt0,Lk+k−1t0,where
(18)L1=ẏ1+A11y1+⋯+A1,kLy1,y2,⋯,yk−1,⋮Lk=Ak,k−1yk−1,Lk+1=Ak−1,k−2yk−2+Ak−1,k−1yk−1,⋮Lk+k−1=s2yk−1+s2yk−2+⋯+sk−2y1,F1=B1u+f1+Lu,u̇,⋯,uk−1,fj,ḟj,⋯,fjk−1,⋮Fk+k−1T=Bk−1uk−2−⋯−1k−1Ak,k−1Ak−1,k−1⋯A21B1+fk−1k−2−⋯−1k−1Ak,k−1Ak−1,k−1⋯A21f1.
This functional guarantees that the critical points ofJy are a solution of (11) and vice versa when y,y~ is selected to be in the nondegenerate bilinear form.Step 4.
The solutions of Equation (13) are the critical points of Equation (11). Using the nondegenerate property on DL,RL, the critical points are a solution to Equation (11). The solution can be approximated using a countable linearly independent (bases) system φ1,φ2,… (complete set of functions; base), where yjt=∑ζ=1naζφζ,n is apositive number; j=1,⋯,k, and aζ is determined so that Jy1t,..,ykt is extremised.Step 5.
J can be derived with respect to aζ,∀ζ=1,⋯,n, where nis a suitable positive number to determine the critical points as a function of aζparameters. One can set dJ/daζ=0,∀ζ=1,..,n to obtain a solvable linear algebraic system. Moreover, ∑ζ=1nΜζ,jaζ=Νj; ∀j=1,⋯,n; M=Μζ,jζ=1,j=1n,k is a nonsingular matrix whose coefficients are obtained from dJ/daζ=0, if yt0∈ω0. Thus, the unique solution is guaranteed, and aζ=Μ−1ζ,jΝj. Otherwise, if yt0 is selected as arbitrary, there may be no solution or infinitely many solutions.Output. The approximate solution yi=∑ζ=1naζφζ,i=1,..,k is obtained.


## 5. Illustration Examples

### 5.1. Example: Hessenberg DAEs of Size 3

For the following differential-algebraic system,(19)ẏ11ẏ12+31−10y11y12+0−211y21y22+1001y31y32=10u+f11f12,(20)ẏ21ẏ22+2101y11y12+1010y21y22=0−1u+f21f22,(21)1001y21y22=10u+f31f32,with ut≜−cost∈c2I,R,
(22)f11=−2t5+t3−t2+5t+4+cost,f12=t5+t4+2t2+t+1,f21=t4+4t3+t2+3t+4,f22=6t4+t2+3t+2−cost,f31=t4+2+cost,f32=t5+t2,A11=31−10,A12=0−211,A13=1001,A21=2101,A22=1010,A32=1001,B1=10,B2=0−1,B3=10y1=y11y12,y2=y21y22,y3=y31y32.Step 1.
AsA32A21A13is an invertible matrix, by differentiating the algebraic constraint with respect to t, we have the following:
(23)A32ẏ2=0−B3u̇+A32B2u−ḟ3+A32f2−A32A21y1−A32A22y2=0.
Equation (23) is a new constraint added to Equation (19). The variable y3 does not appear in Equation (23), and one can differentiate Equation (23) with respect to t to obtain the following expression:
(24)−B3u¨+A32B2u̇−f¨3+A32ḟ2−A32A21ẏ1−A32A22ẏ2=0,⇒−B3u¨+A32B2u̇−f¨3+A32ḟ2−A32A21B1u+f1−A11y1−A12y2−A13y3,−A32A22B2u+f2−A21y1−A22y2=0.
Asc≜A32A21A13 is an invertible matrix, using the implicit function theorem, one can obtain the following:
(25)y3=c−1B3u¨−A32B2u̇+A32A21B1+A32A22B2u+f¨3−A32ḟ2+A32A21f1+A32A22f2−A32A21A11+A32A22A21y1−A32A21A12+A32A22A22y2=Ly1,y2+Lu,u̇,u¨,f1,f2,ḟ2,f¨3,
which means that problem is an index 3 linear Hessenberg DAE. The class of consistent initial conditions can be defined as follows:(26)ω0=y1t0,y2t0∣whichsatisfyA32y2t0=B3ut0+f3t0andA32A21y1t0+A32A22y2t0=−B3u̇t0+A32B2ut0−ḟ3t0+A32f2t0.
From the implicit function theorem andy3=Ly1,y2+Lu,u̇,u¨,f1,f2,ḟ2,f¨3, the reduced-type Hessenberg DAEs can be defined as follows:
(27)ẏ1+A11y1+A12y2+A13Ly1,y2=B1u+f1+Lu,u̇,u¨,f1,f2,ḟ2,f¨3,ẏ2+A21y1+A22y2=B2u+f2,A32y2=B2u+f2,A32A21y1+A32A22y2=−B3u̇+A32B2u−ḟ3+A32f2,A32y2t0=B3ut0+f3t0,A32A21y1t0+A32A22y2t0=−B3u̇t0+A32B2ut0−ḟ3t0+A32f2t0.Step 2.
Based on the previous discussion, one can define the suitable bilinear formy1,y2=y1,Ly2 and determine the functional as follows:
(28)Jω0y1,y2=∫t0tf12L1yTL1y+L2yTL2y+L3yTL3y+L4yTL4y+L3yt0TL3yt0+L4yt0TL4yt0−F1TL1y+F2TL2y+F3TL3y+F4TL4y+F3t0TL3yt0F4t0TL4yt0+dt,where L1y=ẏ1+A11y1+A12y2+A13A~, with A~=−c−1A32A21A11y1+A21A12y2+A22A21y1+A22A22y2.(29)L2y=ẏ2+A21y1+A22y2,L3y=A32y2,L4y=A32A21y1+A32A22y2,F1=B1u−c−1B3u¨+c−1A32B2u̇−c−1A32A21B1+c−1A32A22B2u−c−1f¨3−c−1A32ḟ2+1−c−1A32A21f1−c−1A32A22f2,F2=B2u+f2,F3=B3u+f3,F4=−B3u̇+A32B2u−ḟ3+A32f2.
The variational function with the class of consistent initial conditionsω0 is defined as follows:
(30)Jω0y1,y2=∫t0tf12ẏ1+A11y1+A12y2+A13A~Tẏ1+A11y1+A12y2+A13A~+ẏ2+A21y1+A22y2Tẏ2+A21y1+A22y2+A32y2TA32y2+A32A21y1+A32A22y2TA32A21y1+A32A22y2+A32y2t0TA32y2t0+A32A21y1t0+A32A22y2t0TA32A21y1t0+A32A22y2t0F1Tẏ1+A11y1+A12y2+A13A~+F2Tẏ2+A21y1+A22y2+F3TA32y2+F4TA32A21y1+A32A22y2+F3t0TA32y2t0+F4t0TA32A21y1t0+A32A22y2t0dt.Step 3.
One can set the following linear combination set of independent functions with unknown variablesai,bi,ci, anddi as parameterization solutions over the class of continuously differentiable functions. This study is aimed at determining these unknown parameters. Accordingly, the following condition can be assumed:
(31)y11t=∑i=05aiφit,y12t=∑i=05biφit,y21t=∑i=05ciφit,y22t=∑i=05diφit,where φit=ti, and the unknown parameters ai,bi,ci, and di must be determined using the proposed approach.
Therefore,(32)y31y32=Ly11y12y21y22+Lu,u̇,u¨,f1,f2,ḟ2,f¨3=L∑i=05aiφit,∑i=05biφit,∑i=05ciφit,∑i=05diφit,t+Lu,u̇,u¨,f1,f2,ḟ2,f¨3.
In addition, the functionalJω0y1,y2 becomes a function of unknown coefficients a⟶, b⟶, c⟶, and d⟶.Step 4.
AsJa⟶,b⟶,c⟶,d⟶has a quadratic functional form, to determine its extremals (critical points), one can set∂J/∂a⟶=0, ∂J/∂b⟶=0, ∂J/∂c⟶=0, and ∂J/∂d⟶=0, which leads to the linear algebraic equation Μa⟶b⟶c⟶d⟶=Ν. This problem is directly solvable for a⟶b⟶c⟶d⟶=Μ−1Ν, to obtain the approximate solution y11t,y12t,y21t,y22t,y31t,y32t. The matrix coefficient M is obtained from the derivative of J with respect to a⟶, b⟶, c⟶,and d⟶. Table 1 lists the numerical values of the unknown coefficients a⟶, b⟶, c⟶, and d⟶, after solving the obtained linear algebraic system.Table 1
Unknown coefficients of the representation ofy1t,y2t.
Coeff. ofy11 statesa0a1a2a3a4a5Numerical values110000Coeff. ofy12 statesb0b1b2b3b4b5Numerical values011000Coeff. ofy21 statesc0c1c2c3c4c5Numerical values200010Coeff. ofy2 statesd0d1d2d3d4d5Numerical values001001The numerical zeros of the tables denote anynumber≤10−35. The numerical solutions with equality state constraints are plotted in Figure 1 and are compared with the exact solutions:
(33)y1=1+tt2+t,y2=2+t4t5+t2,y3=t3+t−1+t2,foragivenut=−cost.Figure 1
Comparisons between the proposed approach and the exact solutions.


### 5.2. Example : Hessenberg DAEs of Size 4

Consider the following descriptor system:(34)1001000000001001000000000000000000000000000000001001000000000000ẋ+300021−2412101111102−110010−240000000001−1010000000000−21000010010000x=100−1−1110u+f,where
(35)ut=−e−t,f=6t+e−t+2t3+6t4+4t5+202t5+3t4+4t3+3t+1212+t5+t3+2t2+t3t−e−t+3t2+t3+4t4+5t5+10e−t−1+5t4+3t32t+4+e−t+t3−2t4−2t5t5+t4+e−tt,wheret∈0,1,u∈C3I,R,x=x11,x12,x21,x22,x31,x32,x41,x42T, and f=f11,f122,f21,f22,f31,f32,f41,f42T.Step 1.
To achieve Hessenberg DAEs of size 4, one can rewrite system (34) as follows:
(36a)ẋ1+A11x1+A12x2+A13x3+A14x4=B1u+f1,ẋ2+A21x1+A22x2+A23x3=B2u+f2,ẋ3+A32x2+A33x3=B3u+f3,(36b)A43x3=B4u+f4,with
(37)A11=3000,A12=21−24,A13=102−1,A14=1001,A21=1210,A22=1111,A23=0−240,A32=1−101,A33=00−21,A43=1001,B1=10,B2=0−1,B3=−11,B4=10,f1=f11f12,f2=f21f22,f3=f31f32,f4=f41f42,and the A43A32A21A14 invertible matrix.The variablex4 does not appear in the DAE system; thus, one can differentiate the algebraic constraint (36b) with respect to t and substitute the differential constraints (36a) to obtain the following:
(38)A43ẋ3=B4u̇+ḟ4,A43B3u+f3−A32x2−A33x3=B4u̇+ḟ4.Equation (38) is a new constraint added to the system in Equations (36a) and (36b). The variable x4 does not appear in the last constraint; thus, one can differentiate (38) again with respect to t to obtain the following:
(39)−B4u¨+A43B3u̇−s1B2+s2B3u−f¨4+A43ḟ3−s1f2−s2f3+s1A21x1+s1A22+s2A32x2+s1A23+s2A33x3=0,where s1=A43A32,s2=A43A33.Next, we differentiate Equation (39) for t set:
(40)c≜A43A32A21A14−1=01120.As defined previously,x4 is estimated as follows:
(41)x4=cA43B3u¨+cA43f¨3−cB4u⃛−cf⃛4−cs1B2+cs2B3u̇−cs1ḟ2−cs2ḟ3+cs1A21B1+cs1A22B2+cs1A23B3+cs2A32B2+cs2A33B3u+cs1A21f1+cs1A22+cs2A32f2+cs1A23+cs2A33f3−cs1A21A11+cs1A22A21+cs2A32A21x1−cs1A21A12+cs2A32A22+cs2A33A32+cs1A22A22x2−cs1A21A13+cs1A22A23+cs2A32A23+cs1A23A33+cs2A33A33x3=Lx1,x2,x3+Lu,u̇,u¨,u⃛,f1,f2,f3,ḟ2,ḟ3,f¨3,f⃛4,where the class of the consistent initial conditions is derived as follows:
(42)ω0=x1t0,x2t0,x3t0∣A43x3t0=B4ut0+f4t0−A43A32x2t0−A43A33x3t0=B4u̇t0+ḟ4t0−A43B3ut0−A43f3t0s1A21x1t0+s1A22x2t0+s1A23x3t0+s2A32x2t0+s2A33x3t0=B4u¨t0+f¨4t0−A43B3u̇t0−A43ḟ3t0+s1B2ut0+s1f2t0+s2B3ut0+s2f3t0.Step 2.
The variational function with a class of consistent initial conditions is defined as follows:(43)Jω0x=∫t0tf12ẋ1+A11x1+A12x2+A13x3+A14A~1Tẋ1+A11x1+A12x2+A13x3+A14A~1+ẋ2+A21x1+A22x2+A23x3Tẋ2+A21x1+A22x2+A23x3+ẋ3+A32x2+A33x3Tẋ3+A32x2+A33x3+A43x3TA43x3+−A43A32x2−A43A33x3T−A43A32x2−A43A33x3+A~2TA~2+A43x3t0TA43x3t0+−A43A32x2t0−A43A33x3t0T−A43A32x2t0−A43A33x3t0+A~2t0TA~2t0−F1Tẋ1+A11x1+A12x2+A13x3+A14A~1+F2Tẋ2+A21x1+A22x2+A23x3+F3Tẋ3+A32x2+A33x3+F4TA43x3+F5T−A43A32x2−A43A33x3+F6TA~2+F4t0TA43x3t0+F5t0T−A43A32x2t0−A43A33x3t0+F6t0TA~2t0dt,where
(44)A~1t=−cs1A21A11+cs1A22A21+cs2A32A21x1t−cs1A21A12+cs1A22A22+cs2A32A22+cs2A33A32x2t−cs1A21A13+cs1A22A23+cs1A23A33+cs2A32A23+cs2A33A33x3t,A~2=s1A21x1+s1A22x2+s1A23x3+s2A32x2+s2A33x3,F1=B1u+f1−cA43B3u¨−cA43f¨3+cB4u⃛+cf⃛4+cs1B2u̇+cs1ḟ2−cs1A21B1u−cs1A21f1−cs1A22B2u−cs1A22f2−cs1A23B3u−cs1A23f3+cs2B3u̇+cs2ḟ3−cs2A32B2u−cs2A32f2−cs2A33B3u−cs2A33f3,F2=B2u+f2,F3=B3u+f3,F4=B4u+f4F5=B4u̇+ḟ4−A43B3u−A43f3,F6=B4u¨+f¨4−A43B3u̇−A43ḟ3+s1B2u+s1f2+s2B3u+s2f3.Step 3.
One can setx11t=∑i=05aiφit, x12t=∑i=05biφit, x21t=∑i=05ciφit, x22t=∑i=05diφit, x31t=∑i=05niφit, and x32t=∑i=05miφit, with φit=tito obtain
(45)x4t=Lx1,x2,x3,t+Lu,u̇,u¨,u⃛,f1,f2,f3,ḟ2,ḟ3,f¨3,f⃛4=L∑i=05aiφit,∑i=05biφit,∑i=05ciφit,∑i=05diφit,∑i=05niφit,∑i=05miφit,t+Lu,u̇,u¨,u⃛,f1,f2,f3,ḟ2,ḟ3,f¨3,f⃛4.Step 4.
As discussed earlier, the functionalJ is a function of a⟶, b⟶, c⟶, d⟶, n⟶, and m⟶, and as Ja⟶,b⟶,c⟶,d⟶,n⟶,m⟶ is extremised, one can set ∂J/∂a⟶=0, ∂J/∂b⟶=0, ∂J/∂c⟶=0, ∂J/∂d⟶=0, ∂J/∂n⟶=0, and ∂J/∂m⟶=0, which leads to the linear algebraic equation Μa⟶b⟶c⟶d⟶n⟶m⟶=Ν, which is directly solvable for a⟶b⟶c⟶d⟶n⟶m⟶=Μ−1Ν. Hence, the approximate solution x1t,x2t,x3t,x4t is obtained. The matrix coefficient M is obtained from the derivative of J with respect to a⟶, b⟶, c⟶, d⟶, n⟶, and m⟶. Table 2 presents the numerical values of the unknown coefficients a⟶, b⟶, c⟶, d⟶, n⟶, and m⟶.Table 2
Unknown coefficients of the representation ofx1t,x2t,x3t.
Coeff. ofx11 statesa0a1a2a3a4a5Numerical values410001Coeff. ofx12 statesb0b1b2b3b4b5Numerical values101000Coeff. ofx21 statesc0c1c2c3c4c5Numerical values210000Coeff. ofx22 statesd0d1d2d3d4d5Numerical values310100Coeff. ofx31 statesn0n1n2n3n4n5Numerical values000011Coeff. ofx32 statesm0m1m2m3m4m5Numerical values010000The numerical solutions with equality states are plotted in Figure2 and are compared with the exact solution x1=t5+t+4t2+1,x2=2+tt3+t+3,x3t=t4+t5t,x4t=t3t4+4, for a given ut=−e−t∈C3I,R.Figure 2
Comparison between the proposed approach and the exact solutions.


### 5.3. Example : Hessenberg DAEs of Size 2

The following illustration is taken from [28] and represents the control circuit model of DAEs, which are modeled as Hessenberg linear DAEs of size 2. This illustration is solved using the proposed approach, checking the accuracy of this method using the L2-norm because, to our knowledge, no numerical solution for this illustration exists elsewhere. The following circuit system can be considered:
(46)C100C200000000L000ẋ=00000110−111000R0x+000−1Ue.One can choose the state vectorXt=uC1uC2I2I1, where of all these samples are known and as shown in Table 3. For more details about this circuit model, see [28].Table 3SampleRepresented in a circuit modelUeVoltage source denoting the control inputuC1,uC2Voltages of capacitiesC1,C2I1,I2Amperage of the currents flowing overC1,C2LInductorRResistorStep 1.
Based on these samples, the semiexplicit descriptor system can be rewritten as follows:(47)1001000000001000Ẋ=00000110−11100010X+000−1Ue,(48a)x1̇+A11x1+A12x2=B1Ue,(48b)A21x1=B2Ue,where
(49)x1=uC1uC2I2,x2=I1,A11=00000−11−10,A12=−100,A21=−1,0,−1,B1=000,B2=−1,and (A21A12) is an invertible matrix.This system has index 2 Hessenberg DAEs givenUe∈C1I,R. The variable x2 does not explicitly appear in the algebraic constraint (48b). Thus, by deriving the algebraic constraint (4) with respect to t, the following can be obtained:
(50)A21ẋ1=B2Ue,̇⟹A21B1Ue−A21A11x1−A21A12x2=B2Ue,̇⟹x2=A21A12−1A21B1Ue−A21A11x1−B2Uė.The determinant of (A21A12) is 1; therefore, the inverse exists, and the system has an index of 2. Then, x2=x11t−x12t−B2Uė≜Lx11,x12,x13+Lu,u̇, which is solvable ∀t; Uet=−cost∈C1I,R. The class of consistent initial conditions is defined as follows:
(51)ω0=x11t0,x12t0,x13t0∣−x11t0−x13t0=B2Uet0.Then, the reduced-type Hessenberg DAEs are(52)ẋ11−x11t+x12t=sint,ẋ12−x13=0,ẋ13+x11−x12=0,−x11−x13=B2Ue.Step 2.
The solution is reduced to determine the extremum of the following functional:(53)Jω0x=12L1xt,L1xt+12L2xt,L2xt+12L3xt,L3xt+12L4xt,L4xt+12L4xt0,L4xt0−sint,L1xt−B2Ue,L4xt−B2Uet0,L4xt0=∫t0tf12L1xtTL1xt+L2xtTL2xt+L3xtTL3xt+L4xtTL4xt+L4xt0TL4xt0−sintTL1xt+B2UeTL4xt+B2Uet0TL4xt0dt,where
(54)L1xt=ẋ11−x11t+x12t,L2xt=ẋ12−x13,L3xt=ẋ13+x11−x12,L4xt=−x11−x13.Step 3.
The class of polynomials is a dense set on the class of all continuous functions defined on a compact interval with the maximum norm; thus, the solution can be approximated as a linear combination of linearly independent polynomial functions. Therefore, the following basis functions are assumed for simplicity:(55)x11t=∑i=05aiφit,x12t=∑i=05biφit,x13t=∑i=05ciφit,where φit=ti, and the unknown parameters ai, bi, and ci must be determined using the proposed approach. Then, x2t is determined as follows:
(56)x2t=Lx11,x12,x13,t+Lu,u̇=L∑i=05aiφit,∑i=05biφit,∑i=05ciφit,t+Lu,u̇.Step 4.
One can extremize the following functional:(57)Ja⟶,b⟶,c⟶=∫0112∑i=15aiiti−1−∑i=05aiti+∑i=05bitiT∑i=15aiiti−1−∑i=05aiti+∑i=05biti+12∑i=15biiti−1−∑i=05citi∑i=15biiti−1−∑i=05citi+12∑i=15ciiti−1+∑i=05aiti−∑i=05bitiT∑i=15ciiti−1+∑i=05aiti−∑i=05biti+12−∑i=05aiti−∑i=05citiT−∑i=05aiti−∑i=05citi+12−∑i=05ait0i−∑i=05cit0iT−∑i=05ait0i−∑i=05cit0i−sint∑i=15aiiti−1−∑i=05aiti+∑i=05biti−cost−∑i=05aiti−∑i=05citi−cost0−∑i=05ait0i−∑i=05cit0idt.This process results in∂J/∂a⟶=0, ∂J/∂b⟶=0, and ∂J/∂c⟶=0, which leads to a solvable linear algebraic system with the numerical values of unknown coefficients a⟶, b⟶,andc⟶, as shown in Table 4.Table 4
Numerical coefficients of solutionsx11t,x12t,x13t.
Coeff. ofx11 statesa0a1a2a3a4a5Numerical values11.00141.97890.92800.11410.2397Coeff. ofx12 statesb0b1b2b3b4b5Numerical values0−2.0008−0.4879−0.5545−0.1053−0.1479Coeff. ofx21 statesc0c1c2c3c4c5Numerical values−2−1.0015−1.4784−0.9287−0.2371−0.2371This system can be used in practical applications without an exact solution. Therefore, to test the accuracy of the solution from the present method, one can use theL2-norm by substituting the values of ∑i=05aiφit, ∑i=05biφit, and ∑i=05ciφit into x11t, x12t, and x13t to obtain the following L2-norm error:
(58)ẋ11+A11x11+A12x2−B1Ue+A21x2−B2Ue2=∫01ẋ11−x11+x12−sint2+ẋ12−x132+ẋ13+x11−x122+−x11−x13−cost2dt1/2=7.2617×10−5.The numerical solutions with equality constraints are plotted in Figure3.Figure 3
Comparisons between the proposed approach and the exact solution.As can be seen in theL2-norm, the accuracy of the solutions for the DAEs is less than or equal to 10−5. This accuracy test is excellent and provides the confidence of the solution even without knowing the exact numerical solutions. This test illustration is taken from real-life applications as index 2 DAEs. A simple basis function is applied without any prior information on the nature of the solution types, which proves the reliability and efficiency of the proposed approach.


## 5.1. Example: Hessenberg DAEs of Size 3

For the following differential-algebraic system,(19)ẏ11ẏ12+31−10y11y12+0−211y21y22+1001y31y32=10u+f11f12,(20)ẏ21ẏ22+2101y11y12+1010y21y22=0−1u+f21f22,(21)1001y21y22=10u+f31f32,with ut≜−cost∈c2I,R,
(22)f11=−2t5+t3−t2+5t+4+cost,f12=t5+t4+2t2+t+1,f21=t4+4t3+t2+3t+4,f22=6t4+t2+3t+2−cost,f31=t4+2+cost,f32=t5+t2,A11=31−10,A12=0−211,A13=1001,A21=2101,A22=1010,A32=1001,B1=10,B2=0−1,B3=10y1=y11y12,y2=y21y22,y3=y31y32.Step 1.
AsA32A21A13is an invertible matrix, by differentiating the algebraic constraint with respect to t, we have the following:
(23)A32ẏ2=0−B3u̇+A32B2u−ḟ3+A32f2−A32A21y1−A32A22y2=0.
Equation (23) is a new constraint added to Equation (19). The variable y3 does not appear in Equation (23), and one can differentiate Equation (23) with respect to t to obtain the following expression:
(24)−B3u¨+A32B2u̇−f¨3+A32ḟ2−A32A21ẏ1−A32A22ẏ2=0,⇒−B3u¨+A32B2u̇−f¨3+A32ḟ2−A32A21B1u+f1−A11y1−A12y2−A13y3,−A32A22B2u+f2−A21y1−A22y2=0.
Asc≜A32A21A13 is an invertible matrix, using the implicit function theorem, one can obtain the following:
(25)y3=c−1B3u¨−A32B2u̇+A32A21B1+A32A22B2u+f¨3−A32ḟ2+A32A21f1+A32A22f2−A32A21A11+A32A22A21y1−A32A21A12+A32A22A22y2=Ly1,y2+Lu,u̇,u¨,f1,f2,ḟ2,f¨3,
which means that problem is an index 3 linear Hessenberg DAE. The class of consistent initial conditions can be defined as follows:(26)ω0=y1t0,y2t0∣whichsatisfyA32y2t0=B3ut0+f3t0andA32A21y1t0+A32A22y2t0=−B3u̇t0+A32B2ut0−ḟ3t0+A32f2t0.
From the implicit function theorem andy3=Ly1,y2+Lu,u̇,u¨,f1,f2,ḟ2,f¨3, the reduced-type Hessenberg DAEs can be defined as follows:
(27)ẏ1+A11y1+A12y2+A13Ly1,y2=B1u+f1+Lu,u̇,u¨,f1,f2,ḟ2,f¨3,ẏ2+A21y1+A22y2=B2u+f2,A32y2=B2u+f2,A32A21y1+A32A22y2=−B3u̇+A32B2u−ḟ3+A32f2,A32y2t0=B3ut0+f3t0,A32A21y1t0+A32A22y2t0=−B3u̇t0+A32B2ut0−ḟ3t0+A32f2t0.Step 2.
Based on the previous discussion, one can define the suitable bilinear formy1,y2=y1,Ly2 and determine the functional as follows:
(28)Jω0y1,y2=∫t0tf12L1yTL1y+L2yTL2y+L3yTL3y+L4yTL4y+L3yt0TL3yt0+L4yt0TL4yt0−F1TL1y+F2TL2y+F3TL3y+F4TL4y+F3t0TL3yt0F4t0TL4yt0+dt,where L1y=ẏ1+A11y1+A12y2+A13A~, with A~=−c−1A32A21A11y1+A21A12y2+A22A21y1+A22A22y2.(29)L2y=ẏ2+A21y1+A22y2,L3y=A32y2,L4y=A32A21y1+A32A22y2,F1=B1u−c−1B3u¨+c−1A32B2u̇−c−1A32A21B1+c−1A32A22B2u−c−1f¨3−c−1A32ḟ2+1−c−1A32A21f1−c−1A32A22f2,F2=B2u+f2,F3=B3u+f3,F4=−B3u̇+A32B2u−ḟ3+A32f2.
The variational function with the class of consistent initial conditionsω0 is defined as follows:
(30)Jω0y1,y2=∫t0tf12ẏ1+A11y1+A12y2+A13A~Tẏ1+A11y1+A12y2+A13A~+ẏ2+A21y1+A22y2Tẏ2+A21y1+A22y2+A32y2TA32y2+A32A21y1+A32A22y2TA32A21y1+A32A22y2+A32y2t0TA32y2t0+A32A21y1t0+A32A22y2t0TA32A21y1t0+A32A22y2t0F1Tẏ1+A11y1+A12y2+A13A~+F2Tẏ2+A21y1+A22y2+F3TA32y2+F4TA32A21y1+A32A22y2+F3t0TA32y2t0+F4t0TA32A21y1t0+A32A22y2t0dt.Step 3.
One can set the following linear combination set of independent functions with unknown variablesai,bi,ci, anddi as parameterization solutions over the class of continuously differentiable functions. This study is aimed at determining these unknown parameters. Accordingly, the following condition can be assumed:
(31)y11t=∑i=05aiφit,y12t=∑i=05biφit,y21t=∑i=05ciφit,y22t=∑i=05diφit,where φit=ti, and the unknown parameters ai,bi,ci, and di must be determined using the proposed approach.
Therefore,(32)y31y32=Ly11y12y21y22+Lu,u̇,u¨,f1,f2,ḟ2,f¨3=L∑i=05aiφit,∑i=05biφit,∑i=05ciφit,∑i=05diφit,t+Lu,u̇,u¨,f1,f2,ḟ2,f¨3.
In addition, the functionalJω0y1,y2 becomes a function of unknown coefficients a⟶, b⟶, c⟶, and d⟶.Step 4.
AsJa⟶,b⟶,c⟶,d⟶has a quadratic functional form, to determine its extremals (critical points), one can set∂J/∂a⟶=0, ∂J/∂b⟶=0, ∂J/∂c⟶=0, and ∂J/∂d⟶=0, which leads to the linear algebraic equation Μa⟶b⟶c⟶d⟶=Ν. This problem is directly solvable for a⟶b⟶c⟶d⟶=Μ−1Ν, to obtain the approximate solution y11t,y12t,y21t,y22t,y31t,y32t. The matrix coefficient M is obtained from the derivative of J with respect to a⟶, b⟶, c⟶,and d⟶. Table 1 lists the numerical values of the unknown coefficients a⟶, b⟶, c⟶, and d⟶, after solving the obtained linear algebraic system.Table 1
Unknown coefficients of the representation ofy1t,y2t.
Coeff. ofy11 statesa0a1a2a3a4a5Numerical values110000Coeff. ofy12 statesb0b1b2b3b4b5Numerical values011000Coeff. ofy21 statesc0c1c2c3c4c5Numerical values200010Coeff. ofy2 statesd0d1d2d3d4d5Numerical values001001The numerical zeros of the tables denote anynumber≤10−35. The numerical solutions with equality state constraints are plotted in Figure 1 and are compared with the exact solutions:
(33)y1=1+tt2+t,y2=2+t4t5+t2,y3=t3+t−1+t2,foragivenut=−cost.Figure 1
Comparisons between the proposed approach and the exact solutions.


## 5.2. Example : Hessenberg DAEs of Size 4

Consider the following descriptor system:(34)1001000000001001000000000000000000000000000000001001000000000000ẋ+300021−2412101111102−110010−240000000001−1010000000000−21000010010000x=100−1−1110u+f,where
(35)ut=−e−t,f=6t+e−t+2t3+6t4+4t5+202t5+3t4+4t3+3t+1212+t5+t3+2t2+t3t−e−t+3t2+t3+4t4+5t5+10e−t−1+5t4+3t32t+4+e−t+t3−2t4−2t5t5+t4+e−tt,wheret∈0,1,u∈C3I,R,x=x11,x12,x21,x22,x31,x32,x41,x42T, and f=f11,f122,f21,f22,f31,f32,f41,f42T.Step 1.
To achieve Hessenberg DAEs of size 4, one can rewrite system (34) as follows:
(36a)ẋ1+A11x1+A12x2+A13x3+A14x4=B1u+f1,ẋ2+A21x1+A22x2+A23x3=B2u+f2,ẋ3+A32x2+A33x3=B3u+f3,(36b)A43x3=B4u+f4,with
(37)A11=3000,A12=21−24,A13=102−1,A14=1001,A21=1210,A22=1111,A23=0−240,A32=1−101,A33=00−21,A43=1001,B1=10,B2=0−1,B3=−11,B4=10,f1=f11f12,f2=f21f22,f3=f31f32,f4=f41f42,and the A43A32A21A14 invertible matrix.The variablex4 does not appear in the DAE system; thus, one can differentiate the algebraic constraint (36b) with respect to t and substitute the differential constraints (36a) to obtain the following:
(38)A43ẋ3=B4u̇+ḟ4,A43B3u+f3−A32x2−A33x3=B4u̇+ḟ4.Equation (38) is a new constraint added to the system in Equations (36a) and (36b). The variable x4 does not appear in the last constraint; thus, one can differentiate (38) again with respect to t to obtain the following:
(39)−B4u¨+A43B3u̇−s1B2+s2B3u−f¨4+A43ḟ3−s1f2−s2f3+s1A21x1+s1A22+s2A32x2+s1A23+s2A33x3=0,where s1=A43A32,s2=A43A33.Next, we differentiate Equation (39) for t set:
(40)c≜A43A32A21A14−1=01120.As defined previously,x4 is estimated as follows:
(41)x4=cA43B3u¨+cA43f¨3−cB4u⃛−cf⃛4−cs1B2+cs2B3u̇−cs1ḟ2−cs2ḟ3+cs1A21B1+cs1A22B2+cs1A23B3+cs2A32B2+cs2A33B3u+cs1A21f1+cs1A22+cs2A32f2+cs1A23+cs2A33f3−cs1A21A11+cs1A22A21+cs2A32A21x1−cs1A21A12+cs2A32A22+cs2A33A32+cs1A22A22x2−cs1A21A13+cs1A22A23+cs2A32A23+cs1A23A33+cs2A33A33x3=Lx1,x2,x3+Lu,u̇,u¨,u⃛,f1,f2,f3,ḟ2,ḟ3,f¨3,f⃛4,where the class of the consistent initial conditions is derived as follows:
(42)ω0=x1t0,x2t0,x3t0∣A43x3t0=B4ut0+f4t0−A43A32x2t0−A43A33x3t0=B4u̇t0+ḟ4t0−A43B3ut0−A43f3t0s1A21x1t0+s1A22x2t0+s1A23x3t0+s2A32x2t0+s2A33x3t0=B4u¨t0+f¨4t0−A43B3u̇t0−A43ḟ3t0+s1B2ut0+s1f2t0+s2B3ut0+s2f3t0.Step 2.
The variational function with a class of consistent initial conditions is defined as follows:(43)Jω0x=∫t0tf12ẋ1+A11x1+A12x2+A13x3+A14A~1Tẋ1+A11x1+A12x2+A13x3+A14A~1+ẋ2+A21x1+A22x2+A23x3Tẋ2+A21x1+A22x2+A23x3+ẋ3+A32x2+A33x3Tẋ3+A32x2+A33x3+A43x3TA43x3+−A43A32x2−A43A33x3T−A43A32x2−A43A33x3+A~2TA~2+A43x3t0TA43x3t0+−A43A32x2t0−A43A33x3t0T−A43A32x2t0−A43A33x3t0+A~2t0TA~2t0−F1Tẋ1+A11x1+A12x2+A13x3+A14A~1+F2Tẋ2+A21x1+A22x2+A23x3+F3Tẋ3+A32x2+A33x3+F4TA43x3+F5T−A43A32x2−A43A33x3+F6TA~2+F4t0TA43x3t0+F5t0T−A43A32x2t0−A43A33x3t0+F6t0TA~2t0dt,where
(44)A~1t=−cs1A21A11+cs1A22A21+cs2A32A21x1t−cs1A21A12+cs1A22A22+cs2A32A22+cs2A33A32x2t−cs1A21A13+cs1A22A23+cs1A23A33+cs2A32A23+cs2A33A33x3t,A~2=s1A21x1+s1A22x2+s1A23x3+s2A32x2+s2A33x3,F1=B1u+f1−cA43B3u¨−cA43f¨3+cB4u⃛+cf⃛4+cs1B2u̇+cs1ḟ2−cs1A21B1u−cs1A21f1−cs1A22B2u−cs1A22f2−cs1A23B3u−cs1A23f3+cs2B3u̇+cs2ḟ3−cs2A32B2u−cs2A32f2−cs2A33B3u−cs2A33f3,F2=B2u+f2,F3=B3u+f3,F4=B4u+f4F5=B4u̇+ḟ4−A43B3u−A43f3,F6=B4u¨+f¨4−A43B3u̇−A43ḟ3+s1B2u+s1f2+s2B3u+s2f3.Step 3.
One can setx11t=∑i=05aiφit, x12t=∑i=05biφit, x21t=∑i=05ciφit, x22t=∑i=05diφit, x31t=∑i=05niφit, and x32t=∑i=05miφit, with φit=tito obtain
(45)x4t=Lx1,x2,x3,t+Lu,u̇,u¨,u⃛,f1,f2,f3,ḟ2,ḟ3,f¨3,f⃛4=L∑i=05aiφit,∑i=05biφit,∑i=05ciφit,∑i=05diφit,∑i=05niφit,∑i=05miφit,t+Lu,u̇,u¨,u⃛,f1,f2,f3,ḟ2,ḟ3,f¨3,f⃛4.Step 4.
As discussed earlier, the functionalJ is a function of a⟶, b⟶, c⟶, d⟶, n⟶, and m⟶, and as Ja⟶,b⟶,c⟶,d⟶,n⟶,m⟶ is extremised, one can set ∂J/∂a⟶=0, ∂J/∂b⟶=0, ∂J/∂c⟶=0, ∂J/∂d⟶=0, ∂J/∂n⟶=0, and ∂J/∂m⟶=0, which leads to the linear algebraic equation Μa⟶b⟶c⟶d⟶n⟶m⟶=Ν, which is directly solvable for a⟶b⟶c⟶d⟶n⟶m⟶=Μ−1Ν. Hence, the approximate solution x1t,x2t,x3t,x4t is obtained. The matrix coefficient M is obtained from the derivative of J with respect to a⟶, b⟶, c⟶, d⟶, n⟶, and m⟶. Table 2 presents the numerical values of the unknown coefficients a⟶, b⟶, c⟶, d⟶, n⟶, and m⟶.Table 2
Unknown coefficients of the representation ofx1t,x2t,x3t.
Coeff. ofx11 statesa0a1a2a3a4a5Numerical values410001Coeff. ofx12 statesb0b1b2b3b4b5Numerical values101000Coeff. ofx21 statesc0c1c2c3c4c5Numerical values210000Coeff. ofx22 statesd0d1d2d3d4d5Numerical values310100Coeff. ofx31 statesn0n1n2n3n4n5Numerical values000011Coeff. ofx32 statesm0m1m2m3m4m5Numerical values010000The numerical solutions with equality states are plotted in Figure2 and are compared with the exact solution x1=t5+t+4t2+1,x2=2+tt3+t+3,x3t=t4+t5t,x4t=t3t4+4, for a given ut=−e−t∈C3I,R.Figure 2
Comparison between the proposed approach and the exact solutions.


## 5.3. Example : Hessenberg DAEs of Size 2

The following illustration is taken from [28] and represents the control circuit model of DAEs, which are modeled as Hessenberg linear DAEs of size 2. This illustration is solved using the proposed approach, checking the accuracy of this method using the L2-norm because, to our knowledge, no numerical solution for this illustration exists elsewhere. The following circuit system can be considered:
(46)C100C200000000L000ẋ=00000110−111000R0x+000−1Ue.One can choose the state vectorXt=uC1uC2I2I1, where of all these samples are known and as shown in Table 3. For more details about this circuit model, see [28].Table 3SampleRepresented in a circuit modelUeVoltage source denoting the control inputuC1,uC2Voltages of capacitiesC1,C2I1,I2Amperage of the currents flowing overC1,C2LInductorRResistorStep 1.
Based on these samples, the semiexplicit descriptor system can be rewritten as follows:(47)1001000000001000Ẋ=00000110−11100010X+000−1Ue,(48a)x1̇+A11x1+A12x2=B1Ue,(48b)A21x1=B2Ue,where
(49)x1=uC1uC2I2,x2=I1,A11=00000−11−10,A12=−100,A21=−1,0,−1,B1=000,B2=−1,and (A21A12) is an invertible matrix.This system has index 2 Hessenberg DAEs givenUe∈C1I,R. The variable x2 does not explicitly appear in the algebraic constraint (48b). Thus, by deriving the algebraic constraint (4) with respect to t, the following can be obtained:
(50)A21ẋ1=B2Ue,̇⟹A21B1Ue−A21A11x1−A21A12x2=B2Ue,̇⟹x2=A21A12−1A21B1Ue−A21A11x1−B2Uė.The determinant of (A21A12) is 1; therefore, the inverse exists, and the system has an index of 2. Then, x2=x11t−x12t−B2Uė≜Lx11,x12,x13+Lu,u̇, which is solvable ∀t; Uet=−cost∈C1I,R. The class of consistent initial conditions is defined as follows:
(51)ω0=x11t0,x12t0,x13t0∣−x11t0−x13t0=B2Uet0.Then, the reduced-type Hessenberg DAEs are(52)ẋ11−x11t+x12t=sint,ẋ12−x13=0,ẋ13+x11−x12=0,−x11−x13=B2Ue.Step 2.
The solution is reduced to determine the extremum of the following functional:(53)Jω0x=12L1xt,L1xt+12L2xt,L2xt+12L3xt,L3xt+12L4xt,L4xt+12L4xt0,L4xt0−sint,L1xt−B2Ue,L4xt−B2Uet0,L4xt0=∫t0tf12L1xtTL1xt+L2xtTL2xt+L3xtTL3xt+L4xtTL4xt+L4xt0TL4xt0−sintTL1xt+B2UeTL4xt+B2Uet0TL4xt0dt,where
(54)L1xt=ẋ11−x11t+x12t,L2xt=ẋ12−x13,L3xt=ẋ13+x11−x12,L4xt=−x11−x13.Step 3.
The class of polynomials is a dense set on the class of all continuous functions defined on a compact interval with the maximum norm; thus, the solution can be approximated as a linear combination of linearly independent polynomial functions. Therefore, the following basis functions are assumed for simplicity:(55)x11t=∑i=05aiφit,x12t=∑i=05biφit,x13t=∑i=05ciφit,where φit=ti, and the unknown parameters ai, bi, and ci must be determined using the proposed approach. Then, x2t is determined as follows:
(56)x2t=Lx11,x12,x13,t+Lu,u̇=L∑i=05aiφit,∑i=05biφit,∑i=05ciφit,t+Lu,u̇.Step 4.
One can extremize the following functional:(57)Ja⟶,b⟶,c⟶=∫0112∑i=15aiiti−1−∑i=05aiti+∑i=05bitiT∑i=15aiiti−1−∑i=05aiti+∑i=05biti+12∑i=15biiti−1−∑i=05citi∑i=15biiti−1−∑i=05citi+12∑i=15ciiti−1+∑i=05aiti−∑i=05bitiT∑i=15ciiti−1+∑i=05aiti−∑i=05biti+12−∑i=05aiti−∑i=05citiT−∑i=05aiti−∑i=05citi+12−∑i=05ait0i−∑i=05cit0iT−∑i=05ait0i−∑i=05cit0i−sint∑i=15aiiti−1−∑i=05aiti+∑i=05biti−cost−∑i=05aiti−∑i=05citi−cost0−∑i=05ait0i−∑i=05cit0idt.This process results in∂J/∂a⟶=0, ∂J/∂b⟶=0, and ∂J/∂c⟶=0, which leads to a solvable linear algebraic system with the numerical values of unknown coefficients a⟶, b⟶,andc⟶, as shown in Table 4.Table 4
Numerical coefficients of solutionsx11t,x12t,x13t.
Coeff. ofx11 statesa0a1a2a3a4a5Numerical values11.00141.97890.92800.11410.2397Coeff. ofx12 statesb0b1b2b3b4b5Numerical values0−2.0008−0.4879−0.5545−0.1053−0.1479Coeff. ofx21 statesc0c1c2c3c4c5Numerical values−2−1.0015−1.4784−0.9287−0.2371−0.2371This system can be used in practical applications without an exact solution. Therefore, to test the accuracy of the solution from the present method, one can use theL2-norm by substituting the values of ∑i=05aiφit, ∑i=05biφit, and ∑i=05ciφit into x11t, x12t, and x13t to obtain the following L2-norm error:
(58)ẋ11+A11x11+A12x2−B1Ue+A21x2−B2Ue2=∫01ẋ11−x11+x12−sint2+ẋ12−x132+ẋ13+x11−x122+−x11−x13−cost2dt1/2=7.2617×10−5.The numerical solutions with equality constraints are plotted in Figure3.Figure 3
Comparisons between the proposed approach and the exact solution.As can be seen in theL2-norm, the accuracy of the solutions for the DAEs is less than or equal to 10−5. This accuracy test is excellent and provides the confidence of the solution even without knowing the exact numerical solutions. This test illustration is taken from real-life applications as index 2 DAEs. A simple basis function is applied without any prior information on the nature of the solution types, which proves the reliability and efficiency of the proposed approach.


## 6. Conclusions

A reliable and efficient approach of solving a high-index (index) linear Hessenberg control DAEs is presented in this paper. This approach is recommended for solving any linear differential-algebraic systems (Hessenberg DAEs), the necessary condition problem of linear quadratic differential-algebraic optimal control problem, and linear constrained mechanical systems. For nonlinear Hessenberg differential-algebraic systems, the linearization around the operating points is required before using this approach. The good selection of basis function (spline functions-orthonormal base function, etc.) may improve the numerical accuracy. As one can see from the illustrations, the efficient approximate solutions are obtained even for a low number of polynomial basis functions.



---

*Source: 9706255-2021-10-11.xml*