Selected Applications of the Theory of Connections: A Technique for Analytical Constraint Embedding

In this paper, we consider several new applications of the recently introduced mathematical framework of the Theory of Connections (ToC). This framework transforms constrained problems into unconstrained problems by introducing constraint-free variables. Using this transformation, various ordinary differential equations (ODEs), partial differential equations (PDEs) and variational problems can be formulated where the constraints are always satisfied. The resulting equations can then be easily solved by introducing a global basis function set (e.g., Chebyshev, Legendre, etc.) and minimizing a residual at pre-defined collocation points. In this paper, we highlight the utility of ToC by introducing various problems that can be solved using this framework including: (1) analytical linear constraint optimization; (2) the brachistochrone problem; (3) over-constrained differential equations; (4) inequality constraints; and (5) triangular domains.


Introduction
The Theory of Connections (ToC), introduced by D. Mortari [1] (see also [2][3][4]), is a framework that transforms constrained problems into unconstrained problems by deriving "constrained expressions". Using these expressions, the governing physical models are modified and solved in the space without constraints. Below, we present a brief overview. For univariate functions, Ref. [1] introduced the constrained expression using the form, Licensee  y(x) = g(x) + ∑ k = 1 n η k ℎ k (x), (1) where the h k (x) are n assigned linearly independent functions, g(x) is a free function, and the coefficients η k are derived by imposing the n constraints. The resulting constrained expression always satisfies the constraints, as long as g(x) is defined, where the constraints are defined.
It can be seen that the constraints are always satisfied (y(x 1 ) = y 1 and y′ x 2 = y 2 ′ ), as long as g(x 1 ) and g′ x 2 are defined. This is an example that demonstrates Equation (1). The theory is general and can be used for many other constraints. In the paper, we will present other examples (see also [1,3] for more examples).
In general, the η k coefficients are expressed in terms of the independent variables and in terms of g(x) evaluated where the constraints are defined. The Equation (1) is called constrained expression because it represents all possible functions satisfying the n constraints, regardless of the function g(x). Please note that g(x) can be discontinuous, partially defined, and even the Dirac delta function, as long as there are no constraints specified where the delta function is infinite.
The function y(x) in Equation (1) typically is applied to ODEs or some variational problems. In general, in the example of ODEs, we can consider a problem for y = y(x), subject to the constraints, which can include various types of constraints, (e.g., constraints on the function or derivatives, integral constraints, or so on). Here, F is some equation of y(x) and its derivatives. Using a ToC expression (Equation (1)), this equation is transformed to an unconstrained problem where the unknown parameter is defined by a new function g = g(x). This transformation can be written as, F x, g, g′, …, g (n) = 0, (2) where F is equations completely independent of the problem constraints. Moreover, F differs from F and represents modified constitutive relations. The major advantage of Equation (2) is that it can be solved as an unconstrained optimization problem. Past work has focused on expanding the function of g(x) by a linear basis [5][6][7] such that where h(x) is a vector of global functions and ξ is vector of unknown coefficients used in the optimization process. These coefficients are determined by minimizing the residual of Equation (2) at some collocation points [8,9]. In general, one can use locally supported h(x) functions. One can also avoid a linear functional representation of g(x) and use a non-linear representation, for example, neural networks and support vector machines [10].
The Multivariate Theory of Connections [2] extends the original univariate theory [1] to ndimensions and to any-degree boundary constraints. This extension can be summarized by the expression, where x = {x 1 , x 1 , …, x n } is the vector of n orthogonal coordinates, c(x) is a function specifying the boundary constraints, A(c(x)) is any interpolating function satisfying the boundary constraints, and g(x) is the free function. Several examples of constrained expressions can be found in Refs. [1,2].
The major research effort of Theory of Connections has been applied to solve linear [4] and non-linear [3] ODEs. This has been done by expanding the free function, g(x), in terms of a set of basis functions (e.g., orthogonal polynomials, Fourier, etc.). Linear or iterative nonlinear least-squares is then used to solve for the coefficients of this expansion. This approach to solve ODEs/PDEs has many advantages over traditional methods: (1) it consists of a unified framework to solve IVP, BVP, or multi-values problems, (2) it provides an approximated solution expressed via analytical functions that can be used for subsequent algebraic manipulation, (3) the solution accuracy is usually obtained fast for many application problems, (4) the procedure can be numerically robust (very low condition number), and (5) it can solve ODE/PDE subject to a variety of constraint types: absolute, relative, linear, non-linear, and integral. Additionally, this technique has recently been extended to solve 2-dimensional PDEs [10].
The purpose of this paper is to introduce a set of five new applications of Theory of Connections; applications that are not covered in the previous references and where ToC is found effective. The main motivation and a short rational description of the new applications considered in this paper are:

1.
Analytic linear constraint optimization. In this application, the analytical problem of finding the extreme of a quadratic vectorial form subject to linear constraint is obtained using ToC instead of using the classical Lagrange multipliers technique.

2.
Brachistochrone problem. This application is dedicated to one of the most famous problem in calculus of variations: the brachistochrone problem [11]. This problem is solved numerically and at machine error accuracy via the two-point boundary-value problem of the Euler-Lagrange equation associated with the brachistochrone functional integral definition. In particular, the solution is obtained by deriving the differential equation in polar coordinates.

3.
Over-constrained differential equations. This application takes into consideration the problem of solving ODEs subject to more constraints than the degree of the ODE which arise in many areas of science and engineering. For example, in the problem of orbit determination of a satellite the number of measurements exceeds the order of the governing differential equation [12]. Furthermore, multi-purpose optimization [13][14][15] deals specifically with optimizing across multiple objective functions simultaneous where trade-offs (weighting) between two conflicting objectives must be taken into account. In this section, a weight least-squares solution is provided for over-constrained differential equations by assigning relative weights to the constraints and then solving the weighted constrained ODE by ToC. Additionally, another example showing the sequence of continuous solutions of an IVP morphing into an BVP is provided.

4.
Inequality constraints. This application extends the ToC framework to include inequality constraints [16]. This is obtained using a combination of sigmoid functions to keep the constrained expression within the inequality constraints. This allows for the derivation of functions constrained by user-defined bounds than can be asymmetric, continuous, or symmetric discontinuous functions. The main motivation for this constraint type comes from optimal control problems (bounded control inputs).

Triangular domains.
Validated by mathematical proof, Ref. [2] has extended the original univariate theory [1] to the multivariate theory [2]. This extension represents the multivariate formulation of ToC subject to arbitrary-order derivative constraints in rectangular domains. This provides an analytical procedure to obtain constrained expressions in any orthogonal/rectangular space that can be used to transform constrained problems into unconstrained problems. In Ref. [

Analytic Linear Constraints Optimization
This section shows an analytical application of the Theory of Connections. The problem, known as quadratic programming (QP) subject to equality constraints, consists of deriving a closed-form solution of the problem to find the min/max of the quadratic function, F(x), subject to m < n linear constraints, where Q ∈ ℝ n × n and Q = Q T > 0, c ∈ ℝ n , A ∈ ℝ m × n , and b ∈ ℝ m , are all assigned, and D has rank m (all rows are independent). Let us search a solution in the form, The constraint implies, Therefore, Substituting this expression of x in Equation (3) we obtain, The problem of finding the extreme of Equation (3) is now an unconstrained optimization problem. Stationary condition implies In this example, we show how to transform a constrained optimization problem to unconstrained optimization problem and find the unconstrained g(x), which is the solution of Equation (4). This shows a simple application of Theory of Connections. Unfortunately, the matrix A is singular. However, Equation (4) is a consistent linear system. Therefore, the solution can be found in the null space of A. This is computed by the Moore-Penrose inverse matrix (pseudo-inverse). Let A = CΛC T be the spectral decomposition of A, where the eigenvector matrix, C, is an orthogonal matrix because matrix A is symmetric. Therefore, the solution of the problem given in Equation (3) is, where the diagonal matrix Λ* contains the inverse of the nonzero eigenvalues and zero for zero eigenvalues.

Brachistochrone Problem
The brachistochrone problem is one of many problems arising from the calculus of variations. This body of work focuses on finding the extrema of functionals (mappings from functions to real numbers). Associated functions that maximize or minimize the functionals can be found using the Euler-Lagrange equation. In general, this equation produces an ODE which when solved produces the extremal function. The main application of ToC to the calculus of variation is solving the specific ODEs which can be highly complex.
Consider finding the function, f s (x), minimizing the integral, The solution of this equation, f 0 (x), is an extremal solution (where the derivative is zero, which is a necessary condition). This solution, however, is not sufficient as f 0 (x) may be associated with some local minima, maxima, or saddle points. In other word, after solving Equation (6), we can consider f s (x) = f 0 (x), if the condition, is verified. This can be done, for example, by evaluating the second derivative (Hessian) of F x, f(x), f′(x) . Figure 1 shows three different coordinate systems that may be used to solve the brachistochrone problem.
Typically, the brachistochrone problem is written using coordinate system (a) of Figure 1. In this coordinate system, the total travel time of the sliding bead is given by the functional [11], F x, y, y′ = 1 + y ′2 2gy , y = y(x) (7) subject to, y(0) = 0 and y x 1 = y 1 .
Using the steps of the calculus of variations outlined in Equation (6) on the function in Equation (7) yields the differential equation for the function y = y(x), 2yy″ + y ′2 + 1 = 0. (8) However, this coordinate system produces the case that y′(0) = ∞ at the initial point and cannot be solved using ToC. If the axes are oriented as shown in coordinate system (b) of Figure 1, the functional used in the calculus of variations becomes, F x, y, y′ = 1 + y ′2 2gx , y = y(x) which after performing the steps of the calculus of variations leads to the differential equation for the function y = y(x), 2xy″ + y ′3 + y′ = 0.
The problem with the differential equation given in Equation (9) is that the solution may take on two different values of y(x) for a given value of x. Analytically, the issues with the DEs given in Equations (8) and (9) are circumvented by solving the problems parametrically. Unfortunately, the DE for the brachistochrone written in the parametric form either contains not-a-number or infinite values when implemented numerically, or admits a trivial solution. For example, let x(p) and y(p) each be functions of the parameter p and their derivatives with respect to p be given by the notation: dx/dp = x p and d 2 x/dp 2 = x pp . Then, Equation (9) can be rewritten using the chain rule as 2x y pp x pp + y p x p 3 + y p x p = 0, x = x(p), y = y(p) . (10) Unfortunately, Equation (10) will have infinite values anytime x p = 0 or x pp = 0. Therefore, it cannot be implemented numerically. Alternatively, one could write the DE as 2xy pp x p 3 + x pp y p x p 2 + y p 3 = 0, x = x(p), y = y(p) . (11) However, Equation (11) admits the trivial solution x = 0. Attempts have been made to avoid the trivial solution by giving a starting value for the non-linear least-squares that is far from the solution x = 0, but none of these attempts were successful.
Another option is to express the brachistochrone problem in polar coordinates (coordinate system (c) in Figure 1). If the coordinates are chosen correctly, then the problems of the Cartesian coordinate systems can be avoided, and the problem does not need to be parameterized. A downside of this approach is a significantly more complicated DE; however, by using the ToC method, solving this DE is very easy. Using this coordinate system, the energy of an object sliding without friction is given by, 1 2 mv 2 = mgr sin ϑ, r = r(ϑ) . (12) In addition, the differential path length in polar coordinates is given by, ds = r 2 + r ′2 dϑ, (13) where r′ = dr/dϑ. Furthermore, using Equation (12), the velocity, v, can be written as, v = ds dt = 2grcosϑ . (14) Using Equations (13) and (14) one can express the total time to get from the start point to the end point as Let the integrand of Equation (15) be defined as the functional F ϑ, r(ϑ), r′(ϑ) . Then, the calculus of variations can be used to minimize the travel time. This leads to the DE, after some simplification, given in Equation (16).
Using the differential equation given by Equation (16), the ToC method was used to solve the particular case with initial condition (x 0 , y 0 ) = (0, 0) and final condition (x f , y f ) = (1, 5).

Over-Constrained Differential Equations
The seminal paper on ToC [1] presented Equation (1) as a way to incorporate n constraints to form a constrained expression, this equation is repeated below: This process leads to a (n × n) matrix populated by the values of h k (x) evaluated at the constraint conditions. The values of η k can then be solve by simply inverting this matrix. Yet, through this process, it is not required that number of η coefficients equals the number of applied constraints, and thus constraint matrix can also be written for n values of η and m constraints. Which leads to, Yet, through the ToC method, the constrained expression is formed and then the differential equation is solved. This separation allows for the incorporation of more constraints than the order of the Differential Equation.
The concept of this is highlighted in Figure 3 which provides an outline which distinguishes the ToC approach from classical methods in interpolation and least-squares. Part (c) in Figure 3 shows the original formulation of ToC where the number of the η coefficients was equal to the number of constraints incorporated. Part (d) highlights the work of this section which combines this general interpolation method with a weighted least-squares technique for the constraints.

Second Order Differential Equation with Three Point Constraints. A Numerical Example
As a numerical test, let us consider solving the Differential Equation, y″ + 2y′ + y = 0, y = y(x) .
For this example, let us assume that a dynamical system described by Equation (17)  where the M matrix is defined by imposing the constraints and the W matrix is populated by the scaled variances above. The solution leads to a constrained expression of the form, y(x) = g(x) + β 1 (x) y 1 − g 1 + β 2 (x) y 2 − g 2 + β 3 (x) y 3 − g 3 , such that
In this specific problem, a Monte Carlo simulation of 10,000 trials was conducted to determine space that the function y(x) could occupy given the "observation" uncertainty and subject to dynamics governed by the Differential Equation.
Since the highest order constraint is the first derivative, we can safely define h 1 ( and α is a weight parameter such that α ∈ [0, 1], which is synonymous with the transformation from the IVP to the BVP. From this derivation, V becomes and solving for η leads again Equation (18) of the form, y(x) = g(x) + β 1 (x) y 1 − g 1 + β 2 (x) y 2 − g 2 + β 3 (x) y 3 − g 3 , where the β terms are defined by, Again, g(x) can be defined as g = ξ T h(x) and the unknown coefficient vector ξ can be solved with least-squares. Figure 5 shows this transformation "surface" along with the residuals of the Differential Equation for validation of the method. It can be seen that the mean residual overall α values are on the order of 10 −14 with a standard deviation on the same order. This plot shows the solution of the Differential Equation, y(x), continuously morphing from initial point and derivative constraints to initial and final point constraints.

Inequality Constraints
In some cases, it may be desired to find all possible trajectories between two functions, f u (x) and f d (x), that represent upper and lower continuous bounds, respectively. This means that f u (x) > f d (x), ∀ x ∈ [x min , x max ]. All trajectories between these two functions can be approximated by the following expression, where g(x) represents the free function. In the cases that there are also equality constraints to be applied, the function g(x) can be replaced by Equation (1). Additionally, σ(k, z) = 1 1 + e −kz and dσ(k, z) dz = kσ(k, z)(1 − σ(k, z)), is a sigmoid function with smoothing parameter k. The sigmoid function acts as a "switching function" that subtracts off the portion of g(x) that exceeds the bounds. In this equation, k is not infinite, and thus y(x) will experience some overshooting/undershooting when g(x) > f u (x) or g(x) < f d (x) that must be quantified. To analyze this, let us consider a simple case where the function is only constrained from above with f u (x) and therefore only overshooting exists.
For this case, the overshooting phenomenon is shown in Figure 6.
Since the only free part of the equation is g(x), the maximum value of the entire function, y(x), will be when the derivative ( dy(x)/ dg(x) = 0). Performing this derivative, Equation (20) becomes, which can be solved for g(x). The resulting value, g(x) max , is the value of g(x) that will cause the maximum overshoot in y(x). This process yields, where productlog is the product logarithm function (also-called the omega function or the Lambert W function). Substituting this result back into Equation (19), we can obtain the maximum value that y(x) can take.
A similar derivation can be made for the lower bound f d (x) ; the solution is shown in The functions f u (x) and f d (x) replace f u (x) and f d (x) in Equation (19), resulting in, Johnston et al. Page 14 Equation (21) shows that the value of k is the only user-selected parameter that effects the difference between y(x) max and f u (x), because productlog(e −1 ) is a constant. Increasing the value of k will decrease the difference between y(x) max and f u (x), and also the difference between y(x) min and f d (x). Decreasing the value of k will have the opposite effect. Equation (24) was used to generate Figure 7, which contains randomly generated inequality boundaries and functions. This figure shows that all the random functions stay within the randomly generated inequality boundaries, and thus demonstrates how this technique can be used to enforce inequality boundary constraints.

Triangular Domains
This section demonstrates how to construct a constrained expression with a free function, g(x, y), that will always meet the Dirichlet boundary conditions on a triangular domain. This section is broken down into two subsections. The first proves that an affine transformation exists that can transform the unit triangle into any desired triangle: shown in Figure 8. The second subsection shows how to build the constrained expression for the unit triangle, and how to extend it, via the affine transformation given in the first subsection, to any triangular domain.

Affine Transformation from the Unit Triangle to the Generic Triangle
Consider a unit triangle, which is shown on the left plot of Figure 8. This unit triangle has the vertices defined by the vectors, P 1  x y z = p 3 − p 1 , p 2 − p 1 , The following identities prove the mapping provided by Equation (25), The inverse transformation is also of interest, and can be written as, The inverse transformation will always exist if the triangle is non-degenerate. The first column of A is the vector that points from p 3 to p 1 , and the second column is the vector from that points from p 2 to p 1 . For a triangle, these two vectors will never be parallel. The third column of A is defined as the cross-product of the first two columns, and will therefore be orthogonal to the vectors that make up the first two columns. Therefore, the columns of A are linearly independent, and A is invertible. Therefore, the inverse transform always exists.

ToC Surfaces on the Unit Triangle
Equation (28) can be specified for any surface, Z = G(X,Y). This means that we can write the identity, where G 1 = G(0,Y), G 2 = G(X, 1 − X) or G 2 = G(1 − Y,Y), and G 3 = G(X, 0). Therefore, ToC allows us to write the constrained expression representing all surfaces satisfying the boundary conditions on the unit triangle as, where G(X,Y) is a free function. The constrained expression given in Equation (29) can be written explicitly as,

Conclusions
This work demonstrated some of the applications of ToC as a unified framework. The first application showed how to solve QP problems subject to equality constraints without the use of Lagrange multipliers. Using ToC, the equality constraints are embedded into a constrained expression that transforms the problem into an unconstrained optimization problem. The second application showed how to solve one of the most famous calculus of variations problems, the brachistochrone problem, with ToC. The calculus of variations technique was used to derive a highly non-linear differential equation in polar coordinates, to avoid singularities, which was solved via ToC. The third application explored solving overconstrained differential equations, where there are more constraints than the order of the differential equation. Using ToC through weighted least-squares, we transform the problem into unconstrained problem, which can easily be solved. The fourth application showed how to extend ToC to incorporate inequality constraints by using sigmoid functions. In addition, a discussion on avoiding undershoot/overshoot was included. The fifth and final application showed how to write constrained expressions for triangular domains subject to Dirichlet boundary conditions. Future work will investigate how some of these applications can be extended to solve optimal control problems. The QP problem subject to linear equality constraints and solving the brachistochrone problem represent the first research thrust in that direction, as these are some of the classic problems in the field of optimal control. Moreover, future work will investigate how these applications can be combined with the application on inequality constraints to solve optimal control problems subject to inequality constraints.
Furthermore, future work will investigate applications of over-constrained differential equations. One example is applications in orbit determination where measurements exceed the order of the differential equation governing the dynamics of the orbiting body. By using the approach developed for weighted least-squares solutions of differential equations, the dynamics of the body and variance of observational measurements can be combined into a single optimization problem. Another area of future work includes multi-objective optimization problems.
Additionally, future research will look to extend constrained expressions on triangular domains to include constraints on the derivatives normal to the sides of the triangle. The obvious application of this research extension will be in solving differential equations on triangular domains; however, this research extension also has applications in meshing, domain discretization, and visualization.  Monte Carlo test for 10,000 trials. Part (a) represents the solutions of the differential equation over the varying observed "constraints" with part (b) quantifying the accuracy of the solution. Parts (c)-(e) represent the distribution of the solution values compared with the true value of the "constraints". ToC Triangular Surface with g(x, y) = cos 1 2 expx sin(5y) x.