Volume 2, Number 1 (2017) pp 1-8 doi 10.20448/808.2.1.1.8 | Research Articles
In linear structures routing we have the following problem : find the extremal of given functional, i.e 2D or 3D curve , which must consist of a special type elements. The parameters of elements are limited and their number is unknown. At first we must determine number of elements and after this we can find their optimal parameters.
In the case of 2D extremal we shall consider a broken line and parabolic spline. The broken line is used as a longitudinal profile of railways and the parabolic spline is used as a longitudinal profile of roads. Initial problem was solved as multi-stage process using the methods of dynamic and nonlinear programming. On each stage we consider the different mathematical models of unknown extremal line:
1. A broken line with short elements similar to longitudinal profile of ground. This model allow us to find the initial approximation of unknown line using nonlinear programming.
2. The result of first stage give us opportunity to find number of elements using dynamic programming.
3. We use a special algorithm of nonlinear programming for solving initial problem with fixed number of elements and the result of second stage as initial approximation.
Keywords: Functional, Extremal, Objective function, Nonlinear programming, Dynamic programming, Reduced antigradient.
Citation | Valery I. Struchenkov (2017). Nonlinear and Dynamic Programming Methods for Solving the Variational Problems of a Special Type. Scientific Modelling and Research, 2(1): 1-8.
Copyright: This work is licensed under a Creative Commons Attribution 3.0 License
Funding : This study received no specific financial support.
Competing Interests: The author declares that there are no conflicts of interests regarding the publication of this paper.
History : Received: 28 February 2017/ Revised: 7 April 2017/ Accepted: 18 April 2017/Published: 25 April 2017
Publisher: Online Science Publishing
Many of optimization problems from different areas of practice demand to find the 2D or 3D curve, which satisfy some constraints on its parameters and give minimum of predetermined integral index (a criterion of optimality). In other words we must solve a variational problem with a functional of a special type and some features of unknown extremal. For example we may consider optimal line structure routing (designing of plan and longitudinal profile) as a new construction and reconstruction process. The functional is a cost of road construction or reduced expenditures on construction and operation of the road. The main feature of this problem is that the unknown extremal must consist of elements of a certain type whose parameters are limited.
Segments of straight line, parabola, arcs of circles, clothoids etc. can be used for various types of structures. Number of elements is unknown and all elements must satisfy the constraints on first derivatives and curvature. Railways and auto roads are very expensive constructions, so expediency of optimization of their projects is obvious in particular in the conditions of rugged relief and complex geology.
It was proved at 1975- 1976 years ,when we used the first CAD, in which the design of longitudinal profile of new railways made by optimization program .Marked features of the problem hampering the development of algorithms and programs as well as the lack of interest of designers and builders to reduce construction costs, have led to the fact that the complex mathematical models and founded methods of optimization are not practically used. Considered sufficient to use different kinds of heuristic algorithms in the process of interactive design.
Currently, even in the most developed modern systems of automated design (CAD), such as CARD/1 . Bentley Rail Track ; Topomatic  and Kurilko and Chesheva  computer is used to solve related routine tasks, but not as a tool to develop the optimal project. It is known that in the same conditions, having the same information, different specialists offer a variety of design solutions. Consideration of intuitively appointed variants does not guarantees closeness of end result of this process to the optimum variant. At the same time it is known that small changes of the route position can cause significant changes in the cost of construction . Consequently, the problem of developing adequate mathematical models and mathematically correct algorithms routing of new railways and auto roads remains relevant. Optimization is the main way of CAD development. The purpose of this article - to show how to use a combination of nonlinear and dynamic programming methods for solving the initial variational problem.
The problem is as follows: We must find a three-dimensional curve x (s), y (s), z (s), which gives
Here s is current length and S-unknown total length of the extremal, which connects endpoints A and B. F (x (s), y (s), z (s)) are given function with continuous first partial derivatives with respect to all three their arguments. Function F can contain as arguments the partial derivatives of x (s), y (s), z (s. Additionally may be inequalities constraints on the coordinates of the some points of the extremal. For optimal linear structures routing x (s), y (s) - install the plan of the route, and z (s) - its longitudinal profile. In this case z (s) is a single - valued function, and ,where ν (maximum longitudinal slope ) has a value of a few parts per thousand. Therefore the length of the unknown extremal and the length of its projection on the horizontal plane are practically equal. Constraint on the slope as a rule is bilateral. The curvature in plan and longitudinal profile, as well as the coordinates of individual points (for example, when crossing water sources and existing communications) are limited.
The problem has significant differences from the problems considered in the classic calculus of variations :
- the constraints are inequalities;
- the extremal should consist of elements of a special type.
These features do not allow us to use the classic apparatus .
If functions x (s), y (s) in (1) are given and z (s) unknown , then we have a problem of finding a flat curve
z (s) and z (s) is a single-valued function.
Another particular problem is the search of a curve in the plane XOY i.e. functions x (s), y (s ). In this case, the extremal can be a multivalued function. For example, when designing roads in a mountainous area. For line structures routing search z (s) - corresponds to the optimal design of the longitudinal profile for a given variant of the route plan. Search plane curve x (s), y (s) corresponds to the design route plan in circumstances where the longitudinal profile is uniquely determined by the plan or weakly depends on it. An example is the trenching at a given depth, or design of the reconstructed railway route plan .
Next we will consider some two-dimensional problems for railways routing. If we design the longitudinal profile then the unknown extremal z(s) - is a broken line, which connects two end points. If we denote the ground profile H (s) then the problem is as follows. For given H(s) find a broken z (s), which satisfies all the constraints, and gives
where S0- known length of the plan, and function F corresponds to the cost of construction for elementary length.
Realistic models must take into account the design of transverse sections of subgrade, the presence of culverts and other man-made structures, the distribution of earth masses and production methods of excavation, etc. These models are discussed in detail in [8-10] .
We reduce the problem (2) to a nonlinear programming problem, which has many interesting features, regardless of the specific form of the function F. Since the number of elements required broken line is unknown, we assume that nodes of ground profile and nodes of extremal (i.e. the route longitudinal profile) have the same abscissae. Ground profile is always in the form of broken line with irregular pitch, and this assumption allow us to fix the number of elements n (the dimension of the problem) and fix length si of the elements (in the plan). This yields a broken line with the number of elements more than we need, but due to many constraints its deviation from the final Z (s) is small .
Our plan is following:
Futher we will analytically express all constraints on z (s). For this ame we take as variables zi (i = 1, 2, ..., n), i.e. ordinates of all nodes. We divide all constraints into three groups.
All slopes are small, so we consider the lengths of the elements and their projections to be equal
Integrand (2) becomes the sum [9, 10] and we receive the following problem of nonlinear programming with linear constraints:
Find min Ф (x, c) with Ax£b , where x is unknown vector and c is known vector of parameters. Matrix A and vector b are known. They define a linear system of inequalities. Ф (x, c) is the objective function.
When we design the plan of the route , i.e. find the functions x (s), y (s), we may use segments of straight line, parabola, arcs of circles and clothoids as the elements of the extremal. In this case, we know the minimum permissible length of elements, the maximum allowable curvature and constraints on the coordinates of individual points. The extremal may be represented as a broken line , which consist of more then required but known number of elements. If we take as unknown the coordinates of nodes, then we can present all constraints analytically through these variables. As in the design of the longitudinal profile, we obtain the problem of nonlinear programming. But in this case the system of constraints is nonlinear.
When designing the reconstruction of railway tracks is not convenient to take as unknown the coordinates xi, yi of nodes. Instead of this we use shifts ui along the normals at the points of the existing route plan (CD in Figure 1).
We can express all constraints on the extremal through these shifts ui. Thus we again obtain the problem of nonlinear programming.
For the corresponding constraints we have the following formulas:
1.Curvature k of the circle drawn through any three adjacent points must be within predetermined limits i.e.
kmini ≤k(ui-1,ui, ui+1) ≤ kmaxi.
This constraints can be presented in the form:
Figure-1. Discrete representation of the route
The numerator in (3) includes the Cartesian coordinates of three adjacent nodes, and the denominator is the product of the corresponding distances. Determinant is a double square of the triangle with vertices at these nodes. It is a well-known formula for radius of a circle R = abc / (4S), where a, b, c are the lengths of the sides, S is a square of a triangle.
At first we express the Cartesian coordinates of the intersections of the normales with the route plan through the variables ui (i = 1,2, ..., n). We obtain: xi =x0i +uicos αi; yi =y0i +uisin αi
Here x0i, y0i are the coordinates of the starting point on the i-th normal, αi is the angle of this normal with the axis X. Then using these coordinates we calculate the distances between points and curvature.
2. Constraints on the coordinates of the some points, including intersections other communications and watercourses may be presented very simply: cm ≤ um ≤ dm. For a strictly fixed points cm = dm. Here m is the number of the normal which intersect the communication or watercourse. These constraints may be many.
Nonlinear programming problems, which we obtain using the discrete representation of the extremal, have significant features. These features relate primarily to the system of constraints on the unknown extremal.
All constraints does not contain more than three adjacent variables. This means that in the linear case (the design of the longitudinal profile) matrix of constraints has the block character. The blocks corresponding to the three types constraints mentioned above, contain in each row only one, two or three non-zero elements. So we can organise an iterative process at each step of which, for any set of active constraints, descent direction is determined without solving any systems of linear equations at all. One of the such direction is reduced antigradient [11, 12].
If Ak -matrix of the system of the active constraints at the k-th iteration and N(Ak)- its null - space  the above properties of constraints for any active set make it easy to find a basis in N(Ak) . If C is a basis matrix and f is an antigradient of the objective function, then CCTf, named reduced antigradient, is some descent direction. Indeed, its scalar product on antigradient is nonnegative since (CCTf,f)=(CTf,CTf)≥0.
The system of nonlinear constraints, obtained at a discrete representation of route plan, contains in each of the inequalities no more than three adjacent variables. We can linearize the active constraints in each iteration point and obtain a system of linear equalities, which has important features. This allow us to find a basis in null -space corresponding to the tangent plane. The descent direction from the current iteration point advantageously carry out not on the tangent plane, but along the curve at the boundary surface without disturbing the active nonlinear constraints. Tangent to the curve coincides with the direction of descent in the tangent plane. This is the meaning of of the basic displacements method .
For example, some points with numbers from i-th to (i + m)-th belong the curve of the limited radius (Fig. 2). Its initial position is AB. We need to find such changes for variables ui , ui+1 , …, ui+m, at which the radius the radius doesn't change. This means that the circumference is shifted as a whole, and any displacement can be represented as a combination of two offsets: shift along some straight and rotation in a plane around some center. As the straight for shift we can take any normal with number i ≤ k ≤ i + m, and as the center of rotation we can take the intersection of the normal with the route plan. For example, some point M lies on the normal which number is r. Let a is the displacement along the normal r , and φ is an angle of rotation (Fig.2). If a and φ are known, we calculate the increments ∆uk (i ≤ k ≤ i+m) through a, φ, angles of the normals with OX and uk(i ≤ k ≤ i+m). If there is one point on the circle whose coordinate uj is equal to minimal or maximum value, then it should be taken as the center of rotation, a j-th normal as the line for shift. In this case, a = 0, and all shifts are calculated trough the rotation angle.
Figure-2. Example basic displacement (shift + rotation).
If there are two or more such points on the circle of limited radius, then all changes of coordinates Δuk (i ≤ k ≤ i + m) are zero.
The linearly independent basis displacements are analog of a basis at the linear space. After linearization of the basis displacements we obtain the basis in the corresponding null-space.
The features of the constraints allow us to check the opportunity of modification the active set without complex calculations not only in linear but also in nonlinear case .
The result of the first stage is a broken line consisting of short elements. It satisfies all constraints, except for constraints on the length of elements. Therefore, its deviations from the desired extremal are small . This broken line must be transformed into a sequence of elements of a predetermined type with minimal deviation. The length of all elements must be greater than or equal to the minimum, and all other constraints must be met. We can to use previously developed algorithms of dynamic programming for solving emerging problems of broken line approximation [15, 16].
Characteristically, the constraints in this case contribute to reducing the computational time due to the reduction of the number of allowable variants.
When designing the longitudinal profile of railways broken line consisting of short elements, is converted to another broken line, but with elements with a length not less than a predetermined value. Conversion algorithm is discussed in detail in [8, 15].
In the process of designing highway longitudinal profile we transform the broken line using the second degree parabolas or circle arcs with or without direct inserts. Dynamic programming algorithms for solving this task are given in [8, 16, 17].
When designing a route plan, we without particular difficulties transform the broken line into sequence of circular curves with straight inserts.
Only feature is that the extremal in this case is not an unambiguous function. For that reason we use normales instead of vertical lines . For transformation of the broken line into the following sequence: direct, clothoid, circular, clothoid , direct etc we need in a new algorithm, because dynamic programming for this problem was unsuccessfull.
The result of the broken line transformation is only initial approximation for futher optimization. Now we know the quantity of elements, i.e. dimension of our problem and we have a good initial approximation for using nonlinear programming. At this stage we can consider more realistic models of the objective function and take into account relationship of the elements in cuts and fills, which are constructed together. This relationship does not allow us to calculate the objective function separately for each element and to use dynamic programming . So we need in nonlinear programming algorithm. In the designing of railways final project line consist of linear elements, so we can use the algorithm of the first stage after conversion of the old derivatives of the objective function into its derivatives by the ordinates of endpoints of the relevant elements..
If we use the parabolas for roads design the new variables are the slopes at the endpoints of each element instead of ordinates.
If we put the center of Cartesian coordinate system (l, H) at the beginning of the parabolic element, then:
H=al2 +bl +c, where c and b are the ordinate and the slope at the begin of the element , l is the distance from the begin of the element. So we can calculate the ordinate of any point if we know the slopes bj (j = 1,2, ..., n +1) at the endpoints of elements and the ordinate of the starting point of the profile. First we must calculate all parameters aj using bj ,bj+1 and given lenghts of elements Lj
aj= (bj+1 - bj)/(2 Lj). (4)
cj+1=(bj+1 + bj) Lj /2 + cj (5)
Using formulas (4) and (5) we can calculate all ordinates, i.e. old independent variables, old partial derivatives (gradient) of the objective function and recalculate them into new gradient by slopes. The constraints on the slopes are very simply:
bminj ≤ bj ≤ bmaxj
The constraints on radii:
Lj/Rjmin1 ≤ bj+1 - bj≤ Lj/Rjmin2.
Here Rjmin1 <0 and Rjmin2> 0 the minimum allowable radius of convex and concave curve, respectively.
The system of constraints on slopes and their differences is very simple. For that reason we have a nice opportunity to find descent directions at each iteration of the optimization process. A set of active constraints may changes, so we must not only find descent directions but check opportunity to exclude any constraint from this set.
If unite all strings of matrix A, which corresponds to active constraints (at current iteration) into matrix A1, then descent directions ( vector p) belongs to null-space of A1. Indeed, if we move along descent directions, all active constraints must be active. So A1p=0.
If we know some basis in null-space of A1, i.e. vectors ci (i=1,2,...k), where k is a quntity of the active constraints
( strings of A1) then p=Cy. Here C is a basic matrix ( its columns are vectors ci), and y is the vector of coordinates p in this basis.
If there are no active constraints then null-space is empty and C=E, where E is a unit matrix and vector p is the antigradient.
If some variable (the slope bk) is not included in any active constraint, then corresponding pk=fk. Here vector f is the antigradient. It means that we can find the descent direction separatly for left and right side from point k.
If on some piece of profile all slopes have maximal or minimal value, then we can't change any of them without changing the active set of constraints. It means that corresponding null-space is empty. We will name such pieces as zero sections.
For this case we can easily determine which constrain we may exclude from active set. Indeed, we may exclude any constrain if it will not be violate when we will find the descent direction without it. If we cease to consider the slope br to be limiting, then pr becomes equal to fr, and all other components of the vector p remain zero. It follows that the maximum slope br can be left out from active set for this zero section if fr <0, and the minimum slope br can be left out if fr>0.
If on some piece of profile all differences of the adjacent slope (i.e. the curvatures of the elements) have a limit value, then all corresponding slopes must should receive equal increments, otherwise the active set will change or the constraints will be violated. In this case, a dimensionality of null-space is one and all components of a singl basis vector are equal to one. This situations means that we can only rotate such piece of the profile as a whole around any point. All coorditates of vector p ( a descent direction) in this case are equal.
The base matrix C in this case has one column of ones. Thus, all components of the reduced antigradient
g= CCTf in this case are equal to the sum of the corresponding components of the antigradient f. Indeed, CTf is a number, which is equal to this sum, and CCTf is a vector, all components of which are equal to this sum.
If the descent direction p is a projection f on null-space, then the vector n=f-p is a normal to this null- space, so (n,d)=(f-p,d)=0. Here d is any vector of null- space. But in this case all such vectors have equal coordinates and we can assume that all coordinates of vectod d are equal to one. So (f-p,d)=0 means that all coordinates of vector p are equal to the arithmetic mean of the corresponding coordinates of the antigradient f.
For this case we also can easily determine can we exclude the constraint on difference br+1 -br or no. If we exclude this constraint, then the pk+1 will be equal to fright, i.e. the arithmetic mean of the coordinates of the antigradient f right of k and pk respectively will be equal to fleft. If value br+1 -br is maximal, then we can exclude corresponding constraint if fright - fleft <0, and if br+1 -br is minimal, then we must have fright - fleft <0.
If additionally one slope bs on such piece of limit curvature have a limit value, then null-space become empty and any changes of corresponding slopes are impossible. It is again a zero section. In this case, we may exlude the constraint on this slope if the arithmetic mean of the corresponding coordinates of the antigradient f, i.e. fad<0 and bs is maximal or fad>0 and bs is minimal. If it turns out that the constraint on the slope bs cannot be excluded from active set, then we must check can we exclude the constraints on the differences of the slopes or no. This can be done similarly to a piece of limiting curvature which does not contain the limiting slope. But one of the arithmetic mean will replaced on zero. If we checking the constraint on br+1 -br and slope bs has a limit value, then fleft =0 if r>=s , else fright=0.
If there are more than one limit slope in the piece of limiting curvature, then this is degeneracy.
Significant difficulties in finding the descent direction arise when there are active constraints on the ordinates of the individual points of the extremal. In this case, we can't find the descent direction separately on the pieces between which there are elements without limiting slopes and curvature. The matter is that the ordinate of each point depends on the slopes of all preceding elements. In particular, the condition of reaching the end point connects all slopes in one constraint, regardless of the activity of other constraints.
Early were proposed two ways to overcome this difficulties :
Penalty functions allow for the implementation of a simple algorithm, but they create additional gully and significantly jitter the convergence of the optimization process. However, it is possible to avoid constructing and solving systems of linear equations to find the direction of descent.
As already noted, it is easy to construct a basis matrix in null – space, if there are no active constraints on the ordinates of individual points. Suppose that such matrix is constructed and its columns, that is, the basis vectors are
e1,e2,…,eq. Here q is a number of basis vectors. The active constraint on the ordinate gives an additional row to the matrix of active constraints. Accordingly, the coordinates p1, p2,…,pn of the desired descent vector in the null-space must satisfy the additional condition:
α1p1+α2p2+…+ αk pk=0. ( 6 )
Suppose that a point with a bounded ordinate is on an element with number k at a distance αk from its origin.
Then αi (i=1,2,…,k-1)–the lengths of the preceding elements. Further will regard the vector α, whose first k components are αi(i=1,2,…,k) and the remaining components are zero.
Due to the equality 6, the dimensionality of the null - space has decreased by one and we need to correct the basis.
There are some ways to do it. For example, we take any vector ej , for which (α, ej)≠0. This vector will be exclude from basis and all other vectors we recalculate by formula:
The new vectors are orthogonal α and linearly independent. So it is a new basis.
If there are other active conctraints on ordinates, then formula 7 is applied again.
The extremal may be not single-valued function, for example,it is a curve consisting of straight and circular segments or conjugation of straight and circular by clothoida. In this case we use the lengths of elements and the curvature of the circular curves as variables. This problem was solved by using the DFP-optimization method .
Variational problems in which the desired extremal has features similar to those considered in this paper arise not only in the linear structures routing. Therefore, the combined use of dynamic and nonlinear programming seems promising. Specific models and algorithms should be developed taking into account specific features of the tasks to be solved. Quite so was constructed the new CAD system intended for routing railways and auto road.
The Use, "The use of mathematical optimization techniques and computer engineering at the longitudinal profile of the railways," in Proceedings of the All-Union Scientific Research Institute of Transport Engineering, Transport, Moscow, 1977.
 V. I. Struchenkov, "Mathematical models and optimization in line structure routing: Survey and advanced results," International Journal of Communications, Network and System Sciences. Special Issue Models & Algorithms for Applications, vol. 5, pp. 631-637, 2012. View at Google Scholar | View at Publisher
 V. I. Struchenkov, "Piecewise linear approximation of plane curves with constraints in computer-aided design of railway routes," World Journal of Computer Application and Technology, vol. 2, pp. 1-5, 2014. View at Google Scholar
 V. I. Struchenkov, "Per element approximation of plane curves with re-strictions in computer-aided design of road routes," American Journal of Systems and Software, vol. 1, pp. 20-25, 2013. View at Google Scholar
 V. Struchenkov, "Piecewise parabolic approximation of plane curves with restrictions in computer-aided design of road routes. Transaction on machine learning and artifical intelligence," Society Science and Education. United Kingdom, vol. 1, pp. 15-25, 2013.View at Google Scholar