Index

ABSTRACT

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.

DOI: 10.20448/808.2.1.1.8

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

1. INTRODUCTION

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 [1].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 [2]. Bentley Rail Track [3]; Topomatic [4] and Kurilko and Chesheva [5] 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 [6]. 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.

2. MATHEMATICAL MODEL

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 http://onlinesciencepublishing.com/ ,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 [7]:

- 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 [8].

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 [8].

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.

3. FEATURE  EMERGING  PROBLEMS OF  NONLINEAR  PROGRAMMING  AND THEIR SOLUTION

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 [12] the above properties of constraints for any active set make it easy to find  a basis in N(Ak) [13]. 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 [14].

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 [13].

4. DETERMINING A NUMBER OF ELEMENTS AND  INITIAL  APPROXIMATION

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 [14]. 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 [14]. 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.

5. OPTIMIZATION OF THE INITIAL APPROXIMATION

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 [14]. 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 [14]:

  1. Using the penalty functions for violating  such constraints;
  2. Finding  the descent direction by solving the auxiliary systems of linear equations which number is equal to the number of such active constraints.

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 [11].

6. CONCLUSIONS

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.

REFERENCES

[1]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.

[2]CARD/1, "Retrieved: http://www.card-1.com/en/home/," n.d.

[3] Bentley Rail Track, "Retrieved: " n.d.

[4] R. Topomatic, "Retrieved: http:// www.topomatic.ru," n.d.

[5] K. Ju and V. Chesheva, "Geonics ZhELDOR- SAPR," CADmaster, vol. 1, pp. 12-18, 2007.

[6] S. Yousef and M. J. Shahbazi, "Optimum railway alignment. Retrieved from http://www.uic.org/cdrom/2001/wcrr2001/pdf/sp/2_1_1/210.pdf," n.d.

[7] A. V. Ozhegova and R. G. Nasibullin, "Variacionnoe ischislenie: Zadachi, algoritmy, primery. Kazan," p. 230, 2013.

[8] V. I. Struchenkov, Metody optimizacii trass v SAPR linejnyh sooruzhenij. M., Solon – Press, 2014.

[9] V. I. Struchenkov, "Matematicheskie modeli i metody optimizacii v sistemah proektirovanija trass novyh zheleznyh dorog," Informacionnye Tehnologii, vol. 7, pp. 10-13, 2013.

[10] 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

[11] F. Gill, U. Mjurrej, and M. Rajt, "Prakticheskaja optimizacija: Per. s angl.-M: Mir," p. 509, 1985.

[12] M. Aoki, "Vvedenie v metody optimizacii: Per. s angl.- M.: Nauka," p. 344, 1977.

[13] V. I. Struchenkov, "Nonlinear programming algorithms for CAD systems of line structure routing," World Journal of Computer Application & Technology, vol. 2, pp. 114-120, 2014.

[14] V. I. Struchenkov, Metody optimizacii v prikladnyh zadachah: M: Solon- Press, 2009.

[15] 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 

[16] 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 

[17] 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