Viet David Wu
email: david.wu@lip6.fr

Hung Nguyen

Michel Minoux
email: michel.minoux@lip6.fr

Hai Tran
email: hai.tran@energisme.com

Viet Hung Nguyen
email: vhnguyen@isima.fr

Optimal deterministic and robust selection of electricity contracts

Keywords: electricity contract, convex separable optimization, robust optimization, active constraint set

Introduction

Optimally selecting an energy contract is an important issue for many industrial customers which have an electricity contract for each delivery point. Contracts may need to fix the maximum of total energy consumption in kilowatt-hours and the maximum of peak demand in kilowatt over a large period, typically one year.

The electricity contract problem discussed in [START_REF] Chen | A linear programming approach to the electricity contract capacity problem[END_REF] can be formulated as a linear program and solved as such in polynomial time. In the problem discussed in [START_REF] Chen | A linear programming approach to the electricity contract capacity problem[END_REF], there is a fixed cost when the peak demand does not exceed the contract capacity and some additional cost for excess demand.

Combining the optimal energy contract selection problem with a lot sizing problem by considering renewable energy sources has been investigated in [START_REF] Rodoplu | Energy Contract Optimization for the Single Item Lot Sizing Problem in a Flow-Shop Configuration and Multiple Energy Sources[END_REF].

Using a pricing based on the so-called "subscribed power" tends to enforce people to adapt their electricity consumption behaviour in order to match with the subscription [START_REF] Saele | Subscribed power-testing new power based network tariffs stimulating for demand response[END_REF]. This subscribed power is based on billing the cost of the electricity contract, the cost of some penalty of overpowering with respect to the maximum peak demand defined in the contract and some other fixed costs.

In this paper, we investigate an optimal electricity contract selection problem (ECSP) based on subscribed power which is formulated as minimizing a convex separable function subject to total order constraints. The problem is addressed both in the deterministic case and in the case when consumption data are subject to uncertainty. The paper is organized as follows. Section 2 presents the context and proposes a formulation of the optimization problem considered in its deterministic version. Section 3 presents an efficient algorithm based on an active constraint set approach to solve the optimal selection of electricity contracts. Section 4 investigates a robust version of the problem in which uncertainties on consumption are taken into account. Section 5 presents a series of computational results obtained, both for the deterministic case and for the robust version of the problem. Section 6 presents the main conclusions and some perspectives for future research.

2 Modelling the problem

Mathematical formulation of ECSP

In France, ENEDIS is responsible for the management of 95% of the electricity distribution network. The tariffs for using the electricity grid are referred to as TURPE (see enedis.fr).

Over the year, consumption periods are classified into several predefined categories, where each category may have it own consumption behaviour. Under the 5th version of TURPE, denoted TURPE5, K = 5 categories are considered, namely: winter peak hours [START_REF] Ahuja | Solving the Convex Cost Integer Dual Network Flow Problem[END_REF], winter full hours [START_REF] Ahuja | A Fast Scaling Algorithm for Minimizing Separable Convex Functions Subject to Chain Constraints[END_REF], winter off-peak hours (3), summer full hours (4) and summer off-peak hours [START_REF] Best | Minimizing separable convex functions subject to simple chain constraints[END_REF]. Electricity used during peak hours is billed at a higher rate than electricity used during off-peak hours, creating an incentive for homeowners to reduce energy usage during peak hours. Let T be the total number of 10 minute periods (e.g. T = 52560 for a year) and T i the total number of 10 minutes period in category i. Let M = {1, . . . , 12} be the set of all months of the study period, and ∀i ∈ {1, . . . , K}, M i is the subset of M corresponding to category i. For example, for the winter peak hours category with i = 1, M 1 = {1, 2, 12} which corresponds to the three months of winter. Each small 10 minute period t belongs to a category i and a month m ∈ M i , noted as t ∈ T i,m , and for each t there is a required consumption level c t .

Selecting an electricity contract defined by TURPE5 is to choose for each category of consumption i a value of subscribed power x i ≥ 0 (in kW). The problem of selecting a best contract thus consists in determining a K-dimensional nonnegative vector x = (x 1 , . . . , x K ) minimizing the bill due by the consumer. For each category i the corresponding subscribed power x i represents thresholds of consumption, viewed as a commitment of maximum peak demand for each 10 minutes period.

There is an annual cost for the subscription equal to i s i x i (for given s i > 0 in e /kW/year), this cost is an increasing linear function of x. Excess demand of consumption, called overpowering, is allowed but in that case extra penalty costs are due. The overpowering quantity for t is δ t (x i ) = max(0, c t -x i ). For each category i and month m ∈ M i , the penalty cost is p i t∈Ti,m δ t (x i ) 2 (for given p i > 0 in e /kWh), which is a decreasing function of x i . There is also some fixed cost unrelated to subscribed power (e.g. consumption, transports and taxes) which is not necessary to our analysis. The objective function of the problem is the total annual cost F (x) defined as the sum of the annual costs F i (x i ) for the various categories i ∈ {1, ..., K} :

F i (x i ) = s i x i + p i m∈Mi t∈Ti,m δ t (x i ) 2 . ( 1
)
(see Figure 1 for a graphical representation of a function F i (x i )). In order to ensure the stability of the network load, TURPE5 imposes that subscribed powers in high consumption categories (e.g. winter peak hours) must be lower than subscribed powers in low consumption categories (e.g. summer off-peak hours). This results into a series of K-1 constraints of the form ∀i < K, x i ≤ x i+1 , defining a total order on the x i values.

Let c i (resp.: ci ) be the smallest (resp.: greatest) value of consumed power during the time slots in category i. Then we denote by c (resp.: c) the smallest (resp.: greatest) consumed power among the instance. Then the x i variables should meet the bound constraints 0 ≤ c ≤ x i ≤ c. Let C = c -c + 1 be the number of integer values among which the optimal solution is searched. This parameter C together with the number of categories K will be useful for expressing the worstcase time complexity of our algorithm.

The ECSP can then be formulated as the following mathematical program P:

P : min F (x) = K i=1 F i (x i ) (2)
s.t. x i ≤ x i+1 , 1 ≤ i ≤ K-1, (3)
c ≤ x i ≤ c, 1 ≤ i ≤ K, (4)
The constraints (3) are called the order constraints and the constraints (4) are the bound constraints. A solution for P is a K-dimensional vector x = (x 1 , . . . , x K ) T satisfying the constraints (3) and ( 4). The objective function ( 2) is separable and convex as stated in the following Proposition 1.

Proposition 1 Each F i is convex with respect to x i .
Proof For all i ∈ {1, . . . , T }, let g t (x i ) be an univariate convex function such that

g t (x i ) ≥ 0 for all x i ≥ 0. For all x ∈ IR T , let us denotes G(x) = T i=1 g t (x i ) 2 the Euclidean Norm (L 2 ) of     g 1 (x 1 ) . . . g T (x T )     i.e. G(x) = g 1 (x 1 ) . . . g T (x T ) . Let x =     x1 . . . . xT     ≥ 0, x =     x1 . . . . xT     ≥ 0 and λ ∈ [0, 1].
From the convexity of the g t functions and the fact that the norm is a non decreasing function on the set of non negative vectors, it holds:

       g 1 (λx 1 + (1-λ)x 1 ) ≤ λg 1 (x 1 ) + (1-λ)g 1 (x 1 )
. . .

g T (λx T + (1-λ)x T ) ≤ λg T (x T ) + (1-λ)g T (x T )
Therefore:

g 1 (λx 1 + (1-λ)x 1 ) . . . g T (λx 1 + (1-λ)x 1 ) ≤ λg 1 (x T ) + (1-λ)g 1 (x T ) . . . λg T (x T ) + (1-λ)g T (x T )
Now, from the convexity of the Euclidean norm, it follows:

λg 1 (x 1 ) + (1-λ)g 1 (x 1 ) . . . λg T (x T ) + (1-λ)g T (x T ) ≤ λ g 1 (x 1 ) . . . g T (x T ) + (1-λ) g 1 (x 1 ) . . . g T (x T )
In view of this, the following inequality is shown to hold:

g 1 (λx 1 + (1-λ)x 1 ) . . . g T (λx T + (1-λ)x T ) ≤ λ g 1 (x 1 ) . . . g T (x T ) + (1-λ) g 1 (x 1 ) . . . g T (x T )
and from this, we conclude that:

T i=1 g t (λx i + (1-λ)x i ) 2 ≤ λ T i=1 g t (x i ) 2 + (1-λ) T i=1 g t (x i ) 2 ,
which proves the desires convexity property of G.

For each category i the subscription cost is linear in x i and the penalty cost is of the form:

p i G i (x i ) with G i (x i ) = Ti i=1 g t (x i ) 2
, where g t (x i ) = δ t (x i ) is a convex function of x i . Thus, for each category i, F i is convex w.r.t. x i .

State of the art and motivation

The ECSP can be reformulated as a minimum separable convex cost network optimization problem in a graph featuring K nodes and 2K-2 arcs (see Appendix). Thus, in order to solve P, one might well consider applying some of the existing solution procedures, in particular:

the algorithm proposed in [START_REF] Minoux | A polynomial algorithm for minimum quadratic cost flow problems[END_REF][START_REF] Minoux | Solving integer minimum cost flows with separable convex cost objective polynomially[END_REF] which, in our case, would lead to a worst-case complexity O(K 3 logC); -the algorithms proposed in [START_REF] Ahuja | Solving the Convex Cost Integer Dual Network Flow Problem[END_REF] and in [START_REF] Karzanov | Polynomial Methods for Separable Convex Optimization in Unimodular Linear Spaces with Applications[END_REF] which would lead to a worst-case complexity O(K 2 log(K)log(KC)).

For a survey of solution algorithms for convex cost network optimization problems and related nonlinear optimization problems, we refer the reader to [START_REF] Hochbaum | Complexity and algorithms for convex network optimization and other nonlinear problems[END_REF]. Note that all the previous algorithms use various versions of minimum cost flow network algorithms running on various specific types of graphs. Some of these graphs could be of very big size (e.g. the number of the arcs could be KC in the graph considered in [START_REF] Ahuja | Solving the Convex Cost Integer Dual Network Flow Problem[END_REF]) and the construction of these graphs is not taken into account in the time complexity of algorithm. Moreover, the algorithms in [START_REF] Ahuja | Solving the Convex Cost Integer Dual Network Flow Problem[END_REF] and [START_REF] Karzanov | Polynomial Methods for Separable Convex Optimization in Unimodular Linear Spaces with Applications[END_REF] need advanced data structures such as dynamic trees to achieve the time complexity of O(K 2 log(K)log(KC)). Without them, the worst-case complexity would remain

O(K 3 logC).
The ECSP also belongs to the class of convex separable function minimization problems under total order constraints. This class of problem appears in the context of isotonic regression which is a well-known problem in statistics. Given K data c 1 , . . . , c K , the Isotonic Regression Problem (IRP) consists in finding K values x 1 , . . . , x K minimizing the norm p distance between the solution and the data subject to order constraints x 1 ≤ • • • ≤ x K . Formally, it is formulated as follows:

min p i (x i -c i ) p s.t. x 1 ≤ x 2 ≤ . . . ≤ x K .
The generalized IRP is obtained when the objective function is replaced by any convex separation function. One of the most efficient algorithms for solving generalized IRP has been proposed by Ahuja and Orlin [START_REF] Ahuja | A Fast Scaling Algorithm for Minimizing Separable Convex Functions Subject to Chain Constraints[END_REF] under the name of 'Scaling PAV'. Assuming that computing (1) at any x in R K can be done in constant time O(1), its complexity is O(Klog(C)). In Scaling PAV, instead of calculating the exact optimal solution of each F i (x i ) at each iteration, one only updates the intervals to which the optimal solution should belong. The K intervals are all initially set equal to (c, c) and reduced by one half after each iteration. As new order constraints are successively detected and saturated, the number of intervals is progressively decreased as variables are fusioned. After log 2 ( C ) iterations, all the intervals have reduced length less than , and an -optimal solution is found. Since, in each iteration, one works with a at most K intervals, the worst-case time complexity of Scaling PAV is O(Klog(C)).

Obviously, the Scaling PAV algorithm could be applied to the ECSP addressed here. However, in this context of application, the evaluation of the resulting complexity has to be completely reconsidered, because careful examination reveals that the computational effort required to compute each value F i (x i ) involved in the objective function (2) cannot be assumed to be O(1). Indeed, from the expression given in [START_REF] Ahuja | Solving the Convex Cost Integer Dual Network Flow Problem[END_REF], it is readily seen that computing F i (x i ) for any given value of x i , requires O(N + (x i )) arithmetic operations, where N + (x i ) is defined as the number of time steps t in T i such that c t > x i . Thus, for large values of x i (i.e. values close to max(c t , t ∈ T i ), both N + (x i ) and the computation time are reduced. By contrast, for smaller values of x i , N + (x i ) can be of the same order of magnitude as |T i |, the total number of time steps in the definition of F i , thus leading to possibly much bigger computation time. For further discussion and illustration of this dependence of N + (x i ) on x i , please refer to the forthcoming Section 3.2. Clearly, since |T i | (the number of consumption data in the category i) can be much larger than K, the number of decision variables, this parameter turns out to be a key factor to be taken into account in the evaluation of an algorithm for solving P. Thus, in the search for computational efficiency in the solution of ECSP, the dependence of N + (x i ) (a measure of the computational effort) with respect to x i has to be taken into account.

In the present paper, we propose a special-purpose algorithm for solving ECSP which will turn out to be practically more efficient than Scaling PAV, thanks to a proper exploitation of the specific features of the problem pointed out above. The proposed algorithm, which will be called Optim SP, is based on the same scheme of iterative active constraint set detection as Scaling PAV. However, its implementation makes essential use of an extended asymmetric binary search procedure, thanks to which the computational effort required in each iteration to minimize the relevant part of the objective function (2) can be significantly reduced, as compared with usual binary search. It is also worth mentioning here that this possibility of resorting to this extended binary search procedure could not be applied within the Scaling PAV algorithm for the following reason: when running Scaling PAV, at any given step all the intervals must have the same length, and this cannot be the case in case of using the proposed extended binary search.

3 Optim SP: an efficient algorithm based on an active constraint set approach For k ∈ {1, . . . , K-1}, let the k th order constraint denote the constraint x k ≤ x k+1 . For k ∈ {0, . . . , K-1} let P k denote the relaxation of the problem P where the last K -1 -k order constraints are relaxed, i.e. the order constraints in P k are only the first k order constraints:

P k : min F (x) s.t. x i ≤ x i+1 , 1 ≤ i ≤ k, c ≤ x i ≤ c, 1 ≤ i ≤ K.
Note that P 0 is the relaxation of P without the order constraints and P K-1 is P itself. More generally, for B ⊆ {1, . . . , K -1}, let P(B) denote the relaxation of P where all the order constraints are relaxed except those indexed in B:

P(B) : min F (x) s.t. x i ≤ x i+1 , i ∈ B, c ≤ x i ≤ c, 1 ≤ i ≤ K.
In addition, we denote P = (B) the restriction of P(B) where the order constraints in B are set to equality:

P = (B) : min F (x) s.t. x i = x i+1 , i ∈ B, c ≤ x i ≤ c, 1 ≤ i ≤ K.
Note that, if B = ∅, then P = (B) = P 0 .

Algorithm's overview

Optim SP initially solves P 0 and then performs K-1 iterations thus solving successively P 1 , P 2 , . . . , P K-1 . For k = 0, . . . , K -1, let

x k = (x k 1 , . . . , x k K )
T denote the optimal solution of P k found by Optim SP and let B k denote the set of the indices of the active order constraints in P k associated with x k , i.e. the order constraints in P k that x k satisfies at equality. At each iteration k = 1, . . . , K -1, Optim SP sets B ← B k-1 and x ← x k-1 . Then Optim SP iteratively finds the order constraint of greatest index i in P k violated by x, updates B ← B ∪ {i} and finds a new x by solving P = (B). If no such index can be found then Optim SP sets x k ← x and B k ← B and terminates the iteration k. Optim SP can be formally stated as Algorithm 1.

Algorithm 1 Optim SP 1: Set B ← ∅ 2: Solve P = (B) and set x equal to the obtained optimal solution. 3: for k from 1 to K-1 do 4:

Set b ← T rue. 5: repeat 6:

Find the greatest index 1 ≤ i ≤ k such that x i > x i+1 . 7:

Set b ← F alse if such an index does not exist. 8:

if b=True then 9:

Set B ← B ∪ {i} 10:

Solve P = (B) and update x.

11:

end if 12:

until b=False 13: end for 14: return x Algorithm 2 Ajustable binary search The adjustable binary search starts with a lower bound y and an upper bound ȳ. Until yopt is found, the algorithm computes the integer middle point ymid = (y+a * ȳ)/(1 + a) according to a parameter a, and computes objective values lef t cost and right cost for the two neighbours at ymid-and ymid+ , then it updates either y or ȳ.

1: function yopt ← A Search(F I (y), y, ȳ, a) 2: ← 0.1 3: repeat 4: ymid ← (y + a * ȳ)/(1 + a) 5: lef t cost ← F I (ymid -) 6: right cost ← F I (ymid + ) 7: if lef t cost
Please note that the standard binary search procedure corresponds to the choice of a = 1 in A Search. For this value of a, ymid is the middle of the interval (y, ȳ). For a > 1, ymid is no longer the middle of the interval, and we refer to this case as 'asymmetric binary search'. The motivation behind the use of a > 1 stems from the observation that the computational effort O(N + i (x i )) required to compute the value of one (individual aggregated) component of the objective function as defined in ( 1) is strongly dependent on the value of the parameter x i . This dependence is illustrated in Figure 3 Note that when a > 1 does not allow anymore the binary search to continue on the smaller subinterval on the right of ymid, we set the value of a to 1 (see line 9 of Algorithm 2) and resume the last iterations with the standard symmetric binary search. With this switch on the value of a, as our binary search scheme performs at least as well as the symmetric binary search scheme on all iterations except perhaps one, it can be checked that the number of iterations in A Search is at most log 2 (ȳ -y) + 1.

N + i (xi) for (A) N + i (xi) for (B) N + i (xi) for (C)

Solving P = (B)

To perform line 2 in Optim SP which consists in solving P = (∅), i.e. solving P 0 , we call the A Search function for minimizing each F i (x) for i = 1, . . . , K, i.e. x i = A Search(F i (x), c i , ci , a). It easy to see that after execution of line 2, x is an optimal solution of P = (B) with B = ∅.

We now explain the procedure corresponding to line 10 of Optim SP which consists in solving P = (B) with B = ∅. Let us suppose that we are in iteration k (1 ≤ k ≤ K -1) and we are ready to execute line 10 after having added an index i to B in line 9. Let B = B \ {i}, i.e. the value of B before executing line 9. Let x be the optimal solution of P = (B ), i.e. the value of x before the line 10, we have then x i > x i+1 . Let i and ī respectively denote the smallest and the greatest indices of the block of consecutive indices in B which contains i. Let I = {i, i + 1, . . . , ī, ī + 1}. Line 10 is executed by just calling A Search(F I (x), x ī, x i , a) where F I (x) = j∈I F j (x). This function call returns an integer value yopt and the update of x (in line 10 of Algorithm 1) is done as follows: x j ← x j for j / ∈ I and x j ← yopt for j ∈ I.

Correctness of Optim SP

A partition of the index set {1, 2, . . . , K} (of the variables) J = {J 1 , . . . , J k } is called a block partition if its elements are subsets (or blocks) of consecutive indices. Given an x ∈ R K , x is conform to J if x i = x j for all i and j belonging to the same block in J. The correctness of Optim SP is based on well-known results in the literature. In particular, we quote here the following lemma (Lemma 2 in [START_REF] Ahuja | A Fast Scaling Algorithm for Minimizing Separable Convex Functions Subject to Chain Constraints[END_REF] which uses results in [START_REF] Best | Minimizing separable convex functions subject to simple chain constraints[END_REF]) appropriately restated to fit our context.

Lemma 1 [START_REF] Ahuja | A Fast Scaling Algorithm for Minimizing Separable Convex Functions Subject to Chain Constraints[END_REF] If J = {J 1 , . . . , J k } is a block partition of {1, 2, . . . , K} such that there exists a x ∈ R K conform to J which is also feasible for P then x is an optimal solution for P.

Then we can state:

Theorem 1 The solution x output by Optim SP is optimal for P.

Proof As we cannot find a greatest index 1 ≤ i ≤ K -1 such that x i > x i+1 , x is feasible for P. Let B denote the last value of B after termination of Optim SP then B is the set of active constraints that x satisfies at equality. Let us build a block partition J of {1, 2, . . . , K} from B as follows:

-For each of block of consecutive indices I in B with i and ī as respectively the smallest and the greatest indices in I, let us refer to {i, i + 1, . . . , ī, ī + 1} as a block in J.

-For an index 1 ≤ j ≤ K, if neither j nor j -1 belongs to B, then let us refer to {j} as a single block in J.

It is then easy to see that x is a feasible solution of P conform to J. Hence, by Lemma 1, x is optimal for P.

Complexity of Optim SP

Proposition 2 The worst-case time complexity of Optim SP is O(K 2 log(C)) provided that the evaluation of each F i (x) at a specific point x can be done in O(1).

Proof We can see that the worst case occurs when the set B in line 10 of Optim SP contains successively 1, 2, . . . , K -1 elements and Optim SP calls K -1 times A Search(F I (x), x ī, x i , a) in line 10 with the sets I (as defined in Section 3.3) containing successively 2, 3, . . . , K elements. Hence, the total number of times that the functions

F i (x) (i = 1, . . . , k) are minimized in the executions of A Search is K(K+1) 2 -1.
We should also include the minimization of each F i (x) in the K calls of A Search for solving P 0 at line 2 of Optim SP. Hence, overall there are

K(K+1) 2 -1 + K = K(K+3) 2 -1 minimizations of F i (x) in the calls of A Search.
In each minimization of F i (x) by A Search, F i (x) is evaluated at most log 2 (C) + 1 times as pointed out in Section 3.2. Hence, the number of evaluations of the functions F i (x) is O(K 2 log(C)) in the worst case. . Now, if we consider input data as forecasted consumption, actual realisation of electricity consumption can be anticipated to be different from the forecast. Indeed the deviation from predicted consumption can create greater overpowering, the consequence being to increase the electricity bill. In the sequel, we discuss a robust version of the problem of optimal selection of electricity contract and show how it can still be solved using Optim SP algorithm.

Modelling customer demand

Customer demands are often represented as stochastic process based on influencing factors. Those influencing factors are often external condition (temperature, climate...), physical characteristic of dwelling (type, age...), devices and occupants (occupation and behaviour...), prices and subscriptions [START_REF] Barbier | Modélisation de la consommation électrique à partir de grandes masses de données pour la simulation des alternatives énergétiques du futur[END_REF][START_REF] Kavousian | Determinants of residential electricity consumption: Using smart meter data to examine the effect of climate, building characteristics, appliance stock, and occupants' behavior[END_REF][START_REF] Zhou | Data analysis and stochastic modeling of lighting energy use in large office buildings in China[END_REF]. Works on lightning energy in large office building show that consumption can be accurately simulated by considering occupant behaviour and seasonal variations, and can be described using e.g. Poisson and Normal distribution [START_REF] Zhou | Data analysis and stochastic modeling of lighting energy use in large office buildings in China[END_REF]. The consumption model is useful to predict consumption in the short-term (1 day), medium-term (3 days), long-term (7 days) [START_REF] Tatsa | Modélisation et prévision de la consommation horaire d'électricité au Québec: Comparaison de méthodes de séries temporelles[END_REF] and very long-term (1 year) [START_REF] Barbier | Modélisation de la consommation électrique à partir de grandes masses de données pour la simulation des alternatives énergétiques du futur[END_REF]. Many works are based on precise modelling of each device in household [START_REF] Issi | The Determination of Load Profiles and Power Consumptions of Home Appliances[END_REF][START_REF] Pipattanasomporn | Load Profiles of Selected Major Household Appliances and Their Demand Response Opportunities[END_REF][START_REF] Salerno | An Extreme Learning Machine Approach to Effective Energy Disaggregation[END_REF]. The authors of [START_REF] Durand | Analyse de courbes de consommation électrique par chaines de Markov cachées[END_REF] present an approach based on hidden Markov chains for the modelling and statistical analysis of electric consumption curves.

The "Seasonal and Trend decomposition using Loess" (STL) method [START_REF] Cleveland | STL: a seasonal-trend decomposition[END_REF] decomposes the consumption as: raw data(t) = trend(t) + seasonal(t) + remainder(t).

The 'trend' component corresponds to non-stationary long-term evolution of consumption, the 'seasonal' component corresponds to periodical evolution (e.g. period = 1 week), and the remainder (or the residual) variation is what is left over after fitting the model and will be viewed as a realisation of some random process.

In the following, consumption is viewed as a time series in which each sample is collected every 10 minutes. For each category i, we apply the STL method to decompose the time series and the important parameter to characterize uncertainty is taken to be the standard deviation σ i of the remainder. Then for 10 minute period t ∈ {1, . . . , T i } we will typically consider that consumption can deviate from the trend + seasonal forecast by a value represented as an independent realisation of the appropriate probability distribution. In Section 5 below, it will be assumed that the residuals are independent realizations of some known probabilistic distribution (typically a truncated normal distribution) with zero mean and known standard deviation.

Robust problem formulation

We propose here to define a robust version of the problem using a concept of uncertainty set similar to the one introduced in [START_REF] Bertsimas | The Price of Robustness[END_REF]. Following the latter, the uncertainty set, corresponding to ν uncertain parameters µ 1 , . . . , µν , each taking bounded values between 0 and β, and having a very small probability of all simultaneously taking their maximum value β, would typically be the polyhedron defined by the inequalities 0 ≤ µ i ≤ β and ν i=1 µ i ≤ Γ β, where Γ β is the so-called budget of uncertainty. Increasing Γ β increases the range of possible scenarios against which robustness is to be achieved, thus improving the robustness of the solutions obtained but at the expense of increasing the solution cost.

For each category i, let v = (v 1 , . . . , v |Ti| ) be the variable vector of possible deviations, where deviations typically are bounded as -b t ≤ v t ≤ b t . However, in view of the max operator in the expression of the penalty cost, negative deviations will not increase the objective function value. We consider that for a given category i, all b t (t ∈ T i ) are equal and we denote b i this common value. Finally, we consider each deviation is bounded as 0 ≤ v t ≤ b i , where the value b i is chosen to be of the same order of magnitude as the standard deviation σ i (typically b i = r * σ i with r chosen between 2 and 3).

In line with [START_REF] Bertsimas | The Price of Robustness[END_REF], denoting B i the budget of uncertainty for category i, we allow a total deviation at most equal to B i , thus the deviation v has to satisfy the constraint: Ti t=1 v t ≤ B i . As an indication about how to choose B i , consider n realisations of a random variable normally distributed with zero mean and standard deviation σ i . The sum of those n random realisations has a variance σ 2 = nσ 2 i , i.e. its standard deviation is σ = √ nσ i . In line with this remark, taking V i to be an integer value of the same order as √ T i , we define the uncertainty set U i as:

U i = {v|0 ≤ v t ≤ b i , ∀t, Ti t=1 v t ≤ B i = V i b i }.
For defining the robust version of the problem, the overpowering quantity for period t is now δ r t (x i , v t ) = max(0, (c t +v t )-x i ), note that δ r t (x i , v t ) = δ t (x i -v t ). For given x i , the worst-case penalty cost for category i in the robust version of the problem is obtained as the maximum value over the uncertainty set U i of the penalty cost function, which is defined as:

Π r i (x i ) = max v∈Ui m t∈Ti,m (δ r t (x i , v t )) 2 . ( 5
)
The objective function for the robust formulation is

F r (x) = i F r i (x i ) = i s i x i + p i Π r i (x i )
and the robust version of ECSP can be stated as:

P r : min F r (x) s.t. x i ≤ x i+1 , 1 ≤ i ≤ K-1, c ≤ x i ≤ max 1≤i ≤K { ci + b i }, 1 ≤ i ≤ K.
We note that in the above formulation, the upper bound values on the variables x i have been changed to take into account the possible increase of the maximum value of consumption by b i for category i.

Proposition 3 Each Π r i (x i ) is convex in x i .
Proof For any fixed v ∈ U i , we know from Proposition 1 that the function:

m t∈Ti,m (δ r t (x i , v t )) 2
is convex in x i . Now, since the maximum value in ( 5) is to be determined with respect to the (finite) set of extreme points of U i , Π r i is recognized as the pointwise maximum of a finite collection of functions, each of which is convex in x i . This proves the claimed result.

Corollary 1 Each F r i is convex in x i , and hence P r is a convex program.

Using Optim SP to solve the robust version of the problem

Since it appears that the robust version of the optimum contract selection problem has the same structure as the deterministic version, Optim SP algorithm turns out to be readily applicable to the robust version of the optimum contract selection problem. However in each call to A Search function, for a given category i and a given value x i , whenever we have to evaluate the value F r i (x i ), we need now to solve the sub-problem [START_REF] Best | Minimizing separable convex functions subject to simple chain constraints[END_REF].

Proposition 4 For a given x i , any optimal solution to the sub-problem (5) features V i components equal to b i , all the other components being 0.

Proof For a category i, the objective of the sub-problem is to maximize a convex function w.r.t. v, and it is a well-known fact that the maximum is reached at an extreme point of the polyhedron representing the uncertainty set U i . Now, we just have to observe that, since V i has been assumed to be an integer value, each extreme point of the polyhedron U i features exactly V i non zero components equal to their upper bound values b i .

For the category i and given x i , according to proposition 4 the maximization sub-problem ( 5) is solved when the V i uncertain parameters that maximize the objective function are found i.e. we need to find how to best dispatch those V i time slots over the months M i . This optimum dispatch problem, (5) can be cast into a Multiple-Choice Knapsack Problem (MCKP) with integer variables, which according to [START_REF] Zemel | An O(n) Algorithm for the Linear Multiple Choice Knapsack Problem and Related Problems[END_REF] can be solved in worst-case time complexity O(N ) (N = V i * |M i | is the number of items). Let us explain how this MCKP is formulated.

Let the variable zm,n = 1 if the algorithm assigns n ∈ {0, . . . , V i } deviations to the month m ∈ M i , else zm,n = 0. Let T n i,m denote the subset of t ∈ T i,m achieving the n highest values of c t .

Proposition 5 For a given category i and month m, if exactly n deviation parameters v t , t ∈ T i,m , may take on their upper bound value b i , then the worst-case cost value for month m is obtained when v t = b i , for all t ∈ T n i,m . Proof For a given category i and month m, the penalty cost function is maximized by setting maximum deviations for period t with the highest nonzero δ r t values, i.e. exactly n deviation parameters may take on their upper bound value b i , then the maximum cost value is obtained when v t = b i for all t ∈ T n i,m .

According to proposition 5, when n parameters are set equal to their upper bound b i the partial objective cost for month m is:

γm,n = t∈T n i,m (δ r t (x i , b i )) 2 + t∈Ti,m\T n i,m (δ r t (x i , 0)) 2 .
The MCKP to be solved in order to determine the worst-case cost is then: Recall that M i ⊆ M and V i is an integer value of the same order of magnitude as

MKCP: max

√ T i (≤ √ T )
, thus for a given x i determining the optimal value of ( 5 5 Computational experiments

Experiments on the deterministic version of ECSP

This section is devoted to a detailed computational study of the Optim SP algorithm on a series of typical real instances of the ECSP problem. For each instance, the evolution of CPU time as a function of the value chosen for the a parameter is analyzed. The results obtained are also compared against those which would be obtained by using the Scaling PAV algorithm on the same data.

In our computational experiments, we consider 5 data sets referred to as D 1 to D 5 . Each data set is defined by specifying the 52560 values of electricity consumption of the customer under consideration, recorded in every 10-minute period along a full year (2017). This set of values is decomposed into 5 subsets (time categories) T 1 to T 5 , and each T i is in turn decomposed into several subsets (corresponding to the relevant months during which consumption has been observed (these subsets are denoted T i,m for some m in M ). Data set D 1 corresponds to a big factory in the food industry with subscribed power typically in the range 5000-7500 KW; data set D 2 corresponds to a company providing maintenance of transportation equipment with subscribed power typically in the range 1500-2500 KW; data set D 3 corresponds to an industrial bakery with subscribed power typically in the range 1000-1500 KW; D 4 and D 5 correspond to two big hotels with subscribed power typically in the range 500-1500 KW.

D 1 D 2 [c 1 , ..., c K ] [2269,
The main characteristics of each data set are shown in Table 1, namely:

for each time category, the minimum and maximum values of the c t values for t ∈ T i ; -for each time category i, the minimum value x 0 i of F i (x i ), the component of the objective function corresponding to time category i ; it is observed that in all cases, the resulting 5 values do not meet the order constraints imposed in ECSP; -the five components x * 1 , . . . , x * 5 of the optimum solution x * to the ECSP problem, which can be observed to satisfy the order constraints x i ≤ x i+1 (i = 1, . . . , 4). Table 2 displays the computational results obtained with scaling PAV (first line of the table labelled 'sPAV') and Optim SP (in the following lines for values of a ranging from 1 to 8) for the various instances corresponding to the data sets D 1 to D 5 . For each instance, the column labelled '#evalF' provides the number of calls to the A Search procedure; the column labelled 'Total N + ' provides the the evaluation of the resulting total number of arithmetic operations as measured using the N + i (x i ) values; the column labelled 'cpu' indicates the resulting total cpu time (in milliseconds).

All tests reported have been carried out using Python 3.6.8 on an environment of 8-thread quad-core processor with 16 GB RAM and 2.8 GHz CPU running Windows 10 (64 bits). The results shown in Table 2 suggest the following comments:

in all cases but D 5 , the Optim SP algorithm significantly outperforms Scaling PAV in terms of CPU time for a chosen between 4 and 8, and the comparison remains favorable in the wider range 2 ≤ a ≤ 8. For the data set D 5 , conclusions are less clear, but one can observe that there are still two values of a for which Optim SP outperforms Scaling PAV. -in most cases, it can be observed that there is a fairly wide range of values of a for which Optim SP achieves improved efficiency as compared with one of the best previously known algorithms for solving problems of minimizing a convex separable function under order constraints. This shows that Optim SP features a good deal of robustness with respect to the a parameter, and that the choice of a particular value for the parameter a is not a critical issue: a practical consequence of this is that a good value for a will be easily obtained on an experimental basis, by observing the behavior of the Optim SP algorithm on a sample of a few typical instances of the problem for a few values of a in the range [START_REF] Ahuja | A Fast Scaling Algorithm for Minimizing Separable Convex Functions Subject to Chain Constraints[END_REF][START_REF] Durand | Analyse de courbes de consommation électrique par chaines de Markov cachées[END_REF]. -if, for each result shown in Table 2, one computes the ratio cpu/TotalN + , it is observed that all values obtained typically range between 5e -4 and 7e -4, which shows that cpu is, up to small fluctuations, proportional to TotalN + . This confirms the relevance of our analysis based on the use of TotalN + as a measure of computational effort. Statistical tests on remainders reveal that they are uncorrelated, have zero means and are normally distributed, then by considering a truncated normal distribution between -3σ i and 3σ i will give correct representation of deviation of the consumption. In our experiments, the robust parameter b i is chosen to be 3 * σ i (r = 3), and we will use parameter corresponds to

B i = V i * b i (with V i = q * √ T i )
for different values of q ∈ [0, 0.2, 0.5, 1]. Notice that for q = 0 the deterministic version of the problem is obtained. For the summer off-peak hours, the different values of q considered correspond to the values 0, 26, 65 and 130 for V i .

Numerical results concerning the robust version for the data set D 1

Table 4 provides comparative information about the costs of various optimal or sub-optimal solutions to both the deterministic version (q = 0) and the robust version (q > 0) of ECSP for data set D 1 . Table 4 decomposes into three parts:

the first part of the table concerns the contract currently used by the company (a big factory in the food industry), the corresponding solution (6200, 7000, 7000, 7300, 7300), is called the reference solution; the deterministic cost of this reference solution, as well as its costs in terms of the robust objective function for q = 0.2, 0.5, 1 are provided. -the second part concerns the optimal solutions to the deterministic ECSP (for q = 0) and the robust versions of ECSP (for q = 0.2, 0.5, 1); -the third part provides:

(a) the deterministic cost 111296 of the optimal robust solution obtained for q = 0.5 and q = 1, namely: (5881, 6837, 6837, 7341, 7341); (b) for q = 0.2, 0.5, 1, the robust objective function values of the optimal deterministic solution (5536, 6479, 6479, 6914, 6914). The results in Table 4 (i) the comparison of the costs of the optimal solutions (part 2 of the table) with the costs of the reference solution shows that the latter is suboptimal in all cases: the difference is 4.2% for q = 0, 0.4% for q=0.2, and 0.5% for q = 0.5 and q = 1. It is thus observed that the reference solution turns out to be a fairly good approximation (to within 0.5% or so) of the optimal robust solutions for all the (strictly positive) values of q considered. (ii) the comparison between the last three lines of part 2 and the last three lines of part 3 shows the benefit provided by the optimal robust solution over the optimal deterministic solution in the presence of uncertainty: when q = 0.2, the optimal robust function value is 113750, whereas the robust objective function value for the optimal deterministic solution is 118212, showing that the former leads to an improvement in cos 3.8% over the latter. The same comparison for q = 0.5 and q = 1 would lead to 6.3% and 7.6% improvement respectively. (iii) The comparison between the first line of part 2 of the table and the first line of part 3 shows that the difference between the deterministic objective function value of the optimal robust solution for q = 0.5 and q = 1, and the optimal deterministic solution value (111296-107576= 3720) represents an increase by 3.3%. This value can be interpreted as the price of robustness. (iv) the cpu times required for solving the robust versions of ECSP are about 100 times more than for solving the deterministic version, however they typically do not exceed two or three seconds, which is quite acceptable from the point of view of practical applicability.

To conclude, let us point out that, in spite of the fact that the results in Table 4 only concern the instance D 1 , quite similar conclusions would be obtained by analyzing the other data sets D 2 , . . . , D 5 .

Computational efficiency of Optim SP on the Robust ECSP

In addition to the results discussed in Section 5.2.2, we provide in Table 5 results showing the computational efficiency of Optim SP applied to the robust version of ECSP, in a form similar to Table 2. In these experiments, the value q = 0.5 has been taken. The main comments suggested by the results shown in Table 5 are the following:

since the computational effort required to calculate the value of each component of the objective function of the robust problem is significantly bigger than for the deterministic case, it is all but surprising to observe significantly increased cpu values. On average, they can be seen to be multiplied by a factor in the range 50-100. -concerning the influence of varying the value of the a parameter in the range a = 1 to a = 8, a behavior similar to the one observed for the deterministic case can be noticed: for all instances (except D 5 ) Optim SP outperforms Scaling PAV for a wide range of values for a. -the relationship between cpu and the Total N + measure of computational effort, which was almost linear in the deterministic case, is more intricate in the case of the robust ECSP; this is mainly due to the significant overhead induced by the need of repeatedly solving the MCKP discussed in Section 4.

Conclusion and future works

This paper proposes an efficient polynomial-time algorithm based on an active constraint set approach for the electricity contract selection problem (ECSP). This problem is formulated as minimizing a convex separable function subject to total order constraints. The computational results obtained show that the algorithm features improved CPU time as compared with Scaling PAV [START_REF] Ahuja | A Fast Scaling Algorithm for Minimizing Separable Convex Functions Subject to Chain Constraints[END_REF], one of the best existing methods for this problem.

The robust version of the problem has also been investigated. The construction of an uncertainty set representing realistic scenarios of deviation between realized consumption and forecasted values has been proposed, based on the so-called Seasonal and Trend decomposition method (STL). It has been shown that the robust version of ECSP can be efficiently solved using the same algorithm as for the deterministic version. However, taking uncertainty into account requires computing worst-case values of the objective function. The latter problem is reduced to solving a multiple-choice knapsack problem, which is polynomially solvable.

An interesting direction for future research will be to generalize the approach of robustness. In the present paper, it has been assumed that each consumption of a given category i can only deviate by a value b i which is the same for all the time periods in the category. The next step is to have more detailed evaluation of the σ values within each category, the standard deviation may be different from one month to the next (e.g. the consumption of summer off-peak hours in June may be quite different than the one in September). Another interesting subject for future investigation would be to use more accurate stochastic models, possibly capturing statistical dependence from one time instant to the next to better represent uncertainties.

Fig. 1 :

 1 Fig. 1: Objective function for category i: F i (x i )

3. 2

 2 Adjustable binary search for optimizing univariate convex functions F (x) Each step of Algorithm 1 requires the solution of a problem of univariate convex function minimization over a given interval. To achieve this, we propose an extended binary search procedure called A Search, formally stated as Algorithm 2. A Search is an adjustable version of binary search which depends on the input parameter a > 0. Thus, for given univariate convex function F I and an interval [y, ȳ] where I ⊂ {1, . . . , K} and F I (y) = j∈I F j (y), the function A Search returns: yopt = argmin y≤y≤ȳ y integer F I (y).

  on three typical examples of ECSP corresponding to data sets D 1 and D 5 , two of the data sets used in the computational experiments discussed in Section 3 below. More precisely, Figure 3(A) corresponds to dataset D 1 and i = 1; Figure 3(B) corresponds to data set D 1 and i = 5; Figure 3(C) corresponds to data set D 5 and i = 4. histogram (A) histogram(B) histogram (C)

Fig. 3 :

 3 Fig. 3: Histograms and curves showing the dependence of N + i (x i ) on x i for three examples which are part of the data sets D1 and D5 described in Section 5.

4

  Optim SP for a robust version of the problem

  n ≤ V i , Vi n=0 zm,n ≤ 1, ∀m, zm,n ∈ {0, 1}, ∀m∀n.

Proposition 6

 6 ) takes O(|M | * √ T ) worst-case time complexity. Therefore we can state: The worst-case time complexity of Optim SP applied to the robust version of the optimum contract selection problem is O(K 2 M √ T logC).

Table 1 :

 1 Main characteristics of the data sets D 1 , . . . , D 5

		2056, 72, 1769, 1152]	[1430, 1242, 1189, 910, 358]
	[ c1 , ..., c K ] [5536, 6732, 6667, 7539, 7189]	[2419, 2446, 2307, 1822, 1768]
	x 0	[5536, 6573, 6395, 6949, 6738]	[2419, 2425, 2254, 1730, 1691]
	x *	[5536, 6479, 6479, 6914, 6914]	[2317, 2317, 2317, 2317, 2317]
			D 3	D 4
	[c 1 , ..., c K ]	[533, 284, 208, 183, 2]	[480, 193, 321, 375, 315]
	[ c1 , ..., c K ]	[1514, 1596, 1603, 1610, 1588]	[723, 737, 667, 722, 680]
	x 0	[1514, 1568, 1465, 1541, 1481]	[723, 736, 643, 700, 627]
	x *	[1514, 1526, 1526, 1533, 1533]	[711, 711, 711, 711, 711]
			D 5
		[c 1 , ..., c K ]	[735, 545, 664, 572, 665]
		[ c1 , ..., c K ] [1302, 1464, 1818, 2292, 1594]
		x 0	[1302, 1449, 1238, 1716, 1449]
		x *	[1302, 1382, 1382, 1662, 1662]

Table 2 :

 2 Computational results on the deterministic version of ECSP for the data sets D 1 , . . . , D 5

			D 1				D 2
	a	#evalF	Total N +	cpu	#evalF Total N +	cpu
	sPAV	130	113117	44	110	62766	28
	1	192	103498	50	253	72739	34
	2	174	39913	25	256	29724	15
	3	172	18923	10	234	14463	8
	4	167	11286	6	246	12082	8
	5	166	7754	5	252	15993	9
	6	183	5929	4	248	13183	8
	7	168	4801	5	257	11984	7
	8	169	4037	5	255	10989	7
			D 3				D 4
	a		Total N +	cpu	#evalF Total N +	cpu
	sPAV	110	101162	38	90	93442	41
	1	162	90373	44	199	77807	37
	2	149	49375	24	171	34229	17
	3	149	24729	11	169	18330	9
	4	144	11882	6	164	9760	5
	5	134	6426	4	162	6597	4
	6	139	4248	3	159	22971	11
	7	146	3434	3	160	20586	11
	8	147	2977	2	174	19567	11
					D 5
		a		#evalF Total N +	cpu
		sPAV	110	15022	7
		1		167	23838	12
		2		162	24938	12
		3		155	14982	8
		4		161	11746	6
		5		157	9481	5
		6		160	24505	13
		7		167	22706	12
		8		167	21914	11

Table 3 :

 3 Information on remainder of STL decomposition applied the dataset D 1 5.2 Results with Optim SP on the robust version of ECSP This section presents various computational experiments carried out on the version of ECSP.Information on categories and values of the standard deviation (σ i ) of remainder obtained by applying STL decomposition to the data set D 1 are shown in table3(similar results would be observed for the other data sets). For example, the winter peak hours cover three months from December to February and T i = 1488 10-minute periods, and the standard deviation of remainder is σ i = 116.

	Category	Months	c i	ci	T i	√	T i	σ i
	Winter peak hours	Dec-Feb	2269 5536	1488	38	115
	Winter full hours	Nov-Mar	2056 6732	8688	93	124
	Winter off-peak hours	Nov-Mar	72	6667 11562	107	104
	Summer full hours	Apr-Oct	1769 7539 13920	117	144
	Summer off-peak hours	Apr-Oct	1152 7189 16902	130	134
	5.2.1 Analysis of uncertainties through STL decomposition

Table 4 :

 4 suggest the following comments: Comparing the costs of various solutions to the deterministic and robust versions of ECSP for data set D 1

	q	Subscribed Power (SP)	Cost	cpu
	0	6200, 7000, 7000, 7300, 7300 112275
	0.2 6200, 7000, 7000, 7300, 7300 114274
	0.5 6200, 7000, 7000, 7300, 7300 114315
	1	6200, 7000, 7000, 7300, 7300 114315
	0	5536, 6479, 6479, 7316, 7316 107576	25
	0.2 5881, 6837, 6837, 7287, 7287 113750 2235
	0.5 5881, 6837, 6837, 7341, 7341 113771 2404
	1	5881, 6837, 6837, 7341, 7341 113771 2755
	0	5881, 6837, 6837, 7341, 7341 111296
	0.2 5536, 6479, 6479, 6914, 6914 118212
	0.5 5536, 6479, 6479, 6914, 6914 121448
	1	5536, 6479, 6479, 6914, 6914 123161

Table 5 :

 5 3. Computational results on the robust version of ECSP (q = 0.5) for the data sets D 1 , . . . , D 5 .

			D1				D2
	a	#evalF Total N +	cpu	#evalF Total N +	cpu
	sPAV	130	5733875	3815	110	2781354	2363
	1	198	4786268	4058	213	3537792	3271
	2	173	1398527	2046	183	1067315	1803
	3	181	490332	1534	180	461635	1578
	4	173	223949	1401	176	302938	1406
	5	170	130601	1254	176	178908	1362
	6	182	82507	1327	196	324828	1559
	7	175	55507	1267	185	274221	1457
	8	172	50018	1255	183	243549	1472
			D3				D4
	a	#evalF Total N +	cpu	#evalF Total N +	cpu
	sPAV	110	4399366	3024	100	3949268	2834
	1	163	4121953	3512	201	3074664	2971
	2	148	1900657	2237	184	949165	1832
	3	152	570448	1413	174	295506	1387
	4	151	202383	1282	174	129488	1325
	5	143	106005	1099	172	76850	1279
	6	144	64673	1076	175	56555	1290
	7	139	45094	1036	177	916144	1747
	8	145	39596	1065	185	864616	1757
					D5
		a		#evalF Total N +	cpu
		sPAV	110	539235	1115
		1		171	825202	1695
		2		160	171727	1273
		3		156	643125	1472
		4		156	432818	1388
		5		161	355929	1398
		6		166	1192940	1930
		7		163	1078631	1810
		8		167	997721	1773

Acknowledgements This work was supported by Energisme SAS in supplying data and validating the experiments.

Appendix

We state in this appendix the reformulation of P as a separable convex cost network flow problem.

The graph to be considered contains K nodes numbered 1, 2, . . . , K and a set of 2K -1 arcs, decomposed into K 'primary arcs' and K -1 'secondary arcs'.

the set of primary arcs contains the K -1 arcs of the form (i, i + 1) for i = 1, . . . , K -1, together with a 'return arc' (K, 1). The flow value x i on each of these arcs (including x K for the return arc) has to meet the bound constraints 0 ≤ x i ≤ C (and 0 ≤ x K ≤ C) and the corresponding cost function is F i (F K for the return arc); -the set of secondary arcs is composed of K -1 arcs of the form (1, j) for j = 2, . . . , K. The flow value s j on each of these arcs has to meet the bound constraints 0 ≤ s j ≤ C, and the associated cost function is identically 0.

Problem P then reduces to determining a minimum cost circulation on the above graph. We observe that the flow values s j on the secondary arcs play the role of slack variables for the constraints x j-1 ≤ x j for j = 2, . . . , K.


Professor Mikhail F Sukhinin

Oleg Zubelevich

In the memory of

Keywords: 2010 Mathematics Subject Classification. 70G75, 70F25, 70F20, 70H30, 70H03 Brachistochrone, vakonomic mechanics, holonomic systems, nonholonomic systems, Hamilton principle

We consider different generalizations of the Brachistochrone Problem in the context of fundamental concepts of classical mechanics. The correct statement for the Brachistochrone problem for nonholonomic systems is proposed. It is shown that the Brachistochrone problem is closely related to vakonomic mechanics.

Introduction. The Statement of the Problem

The article is organized as follows. Section 4 is independent on other text and contains an auxiliary material with precise definitions and proofs. This section can be dropped by a reader versed in the Calculus of Variations.

Section 3 is also pure mathematical. In this section we justify the alternative statement of the variational problem for the Brachistochrone with nonholonomic constraints.

Other part of the text is less formal and based on Section 4.

The Brachistochrone Problem was stated by Johann Bernoulli in 1696 and became one of classical problems that gave rise to the Calculus of variations.

Since that time the problem was discussed in different aspects numerous times.

We do not even try to touch this long and celebrated history. This article is devoted to comprehension of the Brachistochrone Problem in terms of the modern Lagrangian formalism and to the generalizations which such a comprehension involves.

In the sequel we assume all the functions to be smooth. Recall that a function is smooth in the closed interval [s 1 , s 2 ] if by definition it belongs to C ∞ (s 1 , s 2 ) and all the derivatives are extended to continuous functions in [s 1 , s 2 ].

This assumption is overly strong but we keep it for simplicity of the wording.

1.1. Holonomic Version of the Problem. Assume we are given with a Lagrangian mechanical system with a kinetic energy T and a potential energy V :

L = T -V, T = 1 2 x T t G(x)x t , V = V (x). (1.1)
By the subscript t we denote the derivative in t. Here x = (x 1 , . . . , x m ) T are the local coordinates on a configuration manifold M and G(x) is the matrix of a positively definite quadratic form, G T (x) = G(x), det G(x) = 0, ∀x ∈ M.

All the functions are smooth. Fix two points x 1 , x 2 ∈ M . There are a lot of curves that connect these two points. Let a curve γ be one of them. Assume that this curve is given by the equation x = x(s), x(s 1 ) = x 1 , x(s 2 ) = x 2 .

(1.2)

So that s is a coordinate in γ.

Impose an additional ideal constraint that makes system (1.1) move along the curve γ only. We obtain a one-degree of freedom system with configuration space γ and the generalized coordinate s. Motion of this system along the curve γ is described by the parameter s such that s = s(t) and x = x(s(t)).

(1.3) Substituting this formula in (1.1) we obtain the Lagrangian of this new one-degree of freedom system:

L (s, s t ) = 1 2 x T s (s)G(x(s))x s (s) s 2 t -V (x(s)).
We want to choose the curve γ so that system (1.1) spends minimal time passing along γ from x 1 to x 2 at the energy level h:

1 2 x T s (s)G(x(s))x s (s) s 2 t + V (x(s)) = h. (1.4)
Let us assume that V (x) < h, ∀x ∈ M.

(1.5) Separating variables in (1.4) we see that the time of passing is given by the formula

τ (γ) = s 2 s 1 x T s (s)G(x(s))x s (s) 2(h -V (x(s)) ds.
Therefore we are looking for a stationary point of the functional τ under the boundary conditions (1.2). By other words, the Brachistochrone curve is a geodesic of the Riemann metric

G(x) 2(h -V (x)) . (1.6)
From proposition 7 (see below) it follows that we can equivalently seek for a stationary point of the functional

T(x(•)) = s 2 s 1 x T s (s)G(x(s))x s (s) 2(h -V (x(s))) ds = s 2 s 1 T h -V ds
with the same boundary conditions. By the Hamilton principle it follows that the Brachistochrone curve is a trajectory of a dynamical system with Lagrangian

L(x, x s ) = T h -V . Introduce a function W = - 1 h -V .
Metric (1.6) presents as follows

G(x) 2(h -V (x)) = 1 2 (-W )G.
By the principle of least action in the Moupertuis-Euler-Lagrange-Jacobi form [START_REF] Arnold | Mathematical Methods of Classical Mechanics[END_REF] we conclude that the Brachistochrone is a trajectory of a system with Lagrangian L = T -W at the zero energy level:

T + W = 0.
1.2. Holonomic Constraints. Assume that the Brachistochrone problem is stated for the system with Lagrangian (1.1) and with ideal constraints w(x) = 0, (

where w = (w 1 , . . . w n ) T (x), n < m is a vector of smooth functions in M such that

rang ∂w ∂x (x) = n, ∀x ∈ M.
This statement brings nothing new: equations (1.7) define a smooth submanifold N in M and all constructed above theory works for the corresponding Lagrangian system without constraints on N .

The Brachistochrone Problem with Differential Constraints

2.1. Discussion. Assume that the Lagrangian system (1.1) is equipped with ideal differential constraints

B(x)x t = 0. (2.1)
Here

B(x) =     b 1 1 (x) b 1 2 (x) • • • b 1 m (x) b 2 1 (x) b 2 2 (x) • • • b 2 m (x) . . . . . . . . . . . . b n 1 (x) b n 2 (x) • • • b n m (x)     is a matrix such that rang B(x) = n < m, ∀x ∈ M.
System (2.1) can equivalently be expressed in terms of differential forms

ω i (x) = b i k (x)dx k , T x = n i=1 ker ω i (x) ⊂ T x M.
Namely, x t (t) belongs to an element T x(t) of the differential system {T x } which is determined by equation (2.1), dim T x = m -n. Recall that if system (2.1) can equivalently be presented in the form (1.7) then it is referred to as a holonomic system. Otherwise (2.1) is a nonholonomic system. What the word "equivalently" means and when such a presentation is possible is a content of the Frobenius theorem [START_REF] Sternberg | Lectures on Differential Geometry[END_REF].

As above we substitute equation (1.3) into Lagrangian (1.1) and separating the variables in the energy integral (1.4) we get the functional τ . Then we notice condition (1.5).

But now we must pick only motions that satisfy constraints (2.1). Substituting (1.3) We postpone the discussion of the boundary conditions for a while but note that by the same reason (proposition 7) we can replace τ with T.

We have obtained the problem of stationary points for the functional T defined on a set of curves x = x(s) that obey constraints (2.2). Such a problem is called the Lagrange problem or the problem of vakonomic mechanics.

This problem provoked a lot of confusion and mistakes in classical mechanics. Many researches thought that such a variational problem is equivalent to the Lagrange-d'Alembert equations for the mechanical system with the Lagrangian L and ideal constraints (2.2) (the variable s plays a role of time). If constrains (2.2) are nonholonomic it is not so.

Here is what in this concern Bloch, Baillieul, Crouch and Marsden write in [START_REF] Bloch | Nonholonomic Mechanics and Control (Interdisciplinary Applied Mathematics)[END_REF]:

It is interesting to compare the dynamic nonholonomic equations, that is, the Lagrange-dAlembert equations with the corresponding variational nonholonomic equations. The distinction between these two different systems of equations has a long and distinguished history going back to the review article of Korteweg [START_REF] Korteweg | Uber eine ziemlich verbreitete unrichtige Behandlungsweise eines Problemes der rollenden Bewegung und insbesondere ber kleine rollende Schwingungen um eine Gleichgewichtslage[END_REF] and is discussed in a more modern context in Arnold, Kozlov, and Neishtadt [START_REF] Arnold | Neishtadt: Dynamical Systems III[END_REF]. (For Kozlovs work on vakonomic systems see, e.g., [START_REF] Kozlov | Realization of nonintegrable constraints in classical mechanics[END_REF] and [START_REF] Kozlov | The problem of realizing constraints in dynamics[END_REF])

As Korteweg points out, there were many confusions and mistakes in the literature because people were using the incorrect equations, namely the variational equations, when they should have been using the Lagrange-dAlembert equations; some of these misunderstandings persist, remarkably, to the present day. The upshot of the distinction is that the Lagrange-dAlembert equations are the correct mechanical dynamical equations, while the corresponding variational problem is asking a different question, namely one of optimal control. Perhaps it is surprising, at least at first, that these two procedures give different equations. What, exactly, is the difference in the two procedures? The distinction is one of whether the constraints are imposed before or after taking variations. These two operations do not, in general, commute. With the dynamic Lagrange-dAlembert equations, we impose constraints only on the variations, whereas in the variational problem we impose the constraints on the velocity vectors of the class of allowable curves.

In case of differential constraints (2.2) the situation with boundary conditions is much more complicated than (1.2).

The main question is as follows. Assume that the points x 1 , x 2 ∈ M are connected with a curve x(s). This curve is a stationary point of τ or T on the set of curves that satisfy constraints (2.2) and connect

x 1 , x 2 .
Is the collection of other smooth paths that connect x 1 , x 2 and satisfy (2.2) large enough to reduce the variational problem to the differential equations? Or by other words, is this collection large enough to construct the Lagrange multipliers method? In general the answer is "no". See remark 1 below.

The author does not know whether the situation will be fixed if we demand the constraints to be completely nonholonomic. Such questions seem to be closely related to the Rashevsky-Chow theorem [START_REF] Rashevsky | Any two points of a totally nonholonomic space may be connected by an admissible line[END_REF], [START_REF] Chow | Uber Systeme von linearen partiellen Differentialgleichungen erster Ordnung[END_REF].

To avoid these hard questions we suggest considering the Brachistochrone Problem with another boundary conditions which guarantee the correct employment of the Lagrange multipliers method in the case of nonholonomic constraints.

The boundary Conditions and the Lagrange multipliers

method. Let Σ ⊂ M be a smooth n-dimensional submanifold such that x ∈ Σ =⇒ T x M = T x Σ ⊕ T x .
For briefness of the story sake we proceed with local constructions. Any point x ∈ Σ has a neighborhood U ⊂ M such that there are local coordinates x in U :

x =          y 1 . . . y n z 1 . . . z m-n          = y z , (2.3)
and

U ∩ Σ = {z = z 2 }, B(x)x s = R(x)y s + Q(x)z s ,
where z 2 is a constant vector; R is an n × n matrix and det R(x) = 0, ∀x ∈ U.

Let us impose the boundary conditions:

x(s 1 ) = x 1 = (y T 1 , z T 1 ) T ∈ U, z(s 2 ) = z 2 . (2.4)
We put no restrictions on y(s 2 ).

The geometric sense of these conditions is clear: the left end of the curve x(s) is nailed at x 1 while the right one can slide along the surface Σ :

x(s 1 ) = x 1 , x(s 2 ) ∈ Σ.

(2.5) Observe that boundary conditions (2.4) differ from ones in the vakonomic mechanics [START_REF] Arnold | Neishtadt: Dynamical Systems III[END_REF]. In vakonomic mechanics the ends of the trajectory are fixed and the variations are not supposed to satisfy the equations of constraints. They satisfy the constraints in some asymptotic sense. Theorem 1. Let x(s) be a stationary point of the functional T on the set of functions x(s) that satisfy (2.2), (2.4).

Then there is a smooth function λ(s) = (λ 1 , . . . , λ n )(s) such that x satisfies the equations d ds

∂L * ∂x s - ∂L * ∂x = 0, L * (s, x, x s ) = L(x, x s ) + λ(s)B(x)x s , (2.6
)

and ∂L ∂y s (x(s 2 ), xs (s 2 )) + λ(s 2 )R(x(s 2 )) = 0. (2.7)
This theorem is a direct consequence from theorem 4. Condition (2.7) can be presented in an invariant form:

∂L * ∂x s (s 2 , x(s 2 ), xs (s 2 ))v = 0, ∀v ∈ T x(s 2 ) Σ.
By proposition 1 the stationary point x preserves the "energy":

xT s (s)G(x(s))x s (s) 2(h -V (x(s))) = const.
Show that system (2.6), (2.2) can be presented in the normal form that is

x ss = Ψ(x, x s , λ), λ s = Λ(x, x s , λ).

(2.8) Indeed, system (2.6) takes the form

x T ss + λ s B = α(x, x s , λ), B = B G h -V -1
.

(2.9)

Differentiate (2.2) to have

Bx ss + γ(x, x s ) = 0.

(2.10)

Substituting the second derivatives from (2.9) to (2.10) we obtain

λ s BB T = ψ(x, x s , λ). (2.11)
It is clear det BB T = 0 and we can express λ s from (2.11) and plug it in (2.9). The way we construct (2.8) shows that the vector-function

f (x, x s ) = B(x)x s
is a first integral of (2.8).

Observe also that equations (2.8) are invariant under the substitution

s → -s, λ → -λ.
If we choose y(s 2 ) and z s (s 2 ) then λ(s 2 ) and y s (s 2 ) are defined from equations (2.7), (2.2); and at least for s 2 -s 1 small, we can solve the Cauchy problem for (2.9), (2.11) backwards. Thus the problem is to find y(s 2 ) and z s (s 2 ) so that the boundary conditions x(s 1 ) = x 1 are satisfied.

These observations prompt the following theorem.

Theorem 2. Assume that x 1 is close enough to Σ and s 2 -s 1 > 0 is small enough. Then the problem (2.6), (2.7), (2.5), (2.2) has a unique smooth solution.

2.3. Proof of Theorem 2. Integrating (2.8) twice we get

y 1 = (s 1 -s 2 )y s (s 2 ) + y(s 2 ) + s 1 s 2 ds s s 2 dξ Ψ y (x(ξ), x s (ξ), λ(ξ)), z 1 = (s 1 -s 2 )z s (s 2 ) + z 2 + s 1 s 2 ds s s 2 dξ Ψ z (x(ξ), x s (ξ), λ(ξ)).
(2.12) Here Ψ = (Ψ T y , Ψ T z ) T . Without loss of generality put x 1 = 0. Then the closeness of x 1 to Σ means that |z 2 | is small.

Taking into account (2.2) rewrite (2.12) as

0 = -(s 1 -s 2 )R -1 (x(s 2 ))Q(x(s 2 ))z s (s 2 ) + y(s 2 ) + s 1 s 2 ds s s 2 dξ Ψ y (x(ξ), x s (ξ), λ(ξ)), 0 = (s 1 -s 2 )z s (s 2 ) + z 2 + s 1 s 2 ds s s 2 dξ Ψ z (x(ξ), x s (ξ), λ(ξ)).
(2.13) Introduce notations:

ε = s 2 -s 1 , z 2 = εẑ 2 , y(s 2 ) = εα 1 , z s (s 2 ) = α 2 ;
and put

C(ε, εα 1 ) = -R -1 (x)Q(x) x=ε(α T 1 ,ẑ T 2 ) T .
System (2.13) takes the form

F 1 (ε, α 1 , α 2 ) = -C(ε, εα 1 )α 2 + α 1 + O(ε) = 0, F 2 (ε, α 1 , α 2 ) = -α 2 + ẑ2 + O(ε) = 0. (2.14)
Here we little bit informally write

s 1 s 2 ds s s 2 dξ Ψ(x(ξ), x s (ξ), λ(ξ)) = O(ε 2 ).
The proof is accomplished by application of the Implicit Function Theorem to (2.14):

det ∂F i ∂α j = 0, ε = 0, α 2 = ẑ2 , α 1 = C(0, 0)ẑ 2 .

Minimum of the Functional T

Here we prove that if x 1 is close enough to Σ then the functional T attains a minimum in the class of curves x(s) such that (2.5) and (2.2) holds. The uniqueness follows from theorem 2.

In this section all the inessential positive constants we denote by c, c 1 , c 2 , . . ..

3.1.

The Main Spaces. For details on the Sobolev spaces see [START_REF] Adams | Sobolev Spaces[END_REF].

Let H 1 (s 1 , s 2 ) be the standard Sobolev space of vectors

x(s) = (x 1 , . . . , x m ) T (s) such that

x 2 H 1 (s 1 ,s 2 ) = m i=1 s 2 s 1 (x i (s)) 2 + (x i s (s)) 2 ds < ∞. (3.1)
Recall that the space

H 1 (s 1 , s 2 ) is compactly embedded in C[s 1 , s 2 ].
Introduce a subspace

H 0 = {x ∈ H 1 (s 1 , s 2 ) | x(s 1 ) = 0}.
A function

x 2 = m i=1 s 2 s 1 x i s (s) 2 ds
is a norm in H 0 and this norm is equivalent to the norm (3.1). Moreover,

x C[s 1 ,s 2 ] ≤ c 1 √ s 2 -s 1 x . (3.2) Let a set J ⊂ H 0 consist of functions x(s) = (y T , z T ) T (s) such that z(s 2 ) = z 2
and for almost all s ∈ (s Take a point x ∈ U ∩ Σ such that V (x ) < h. Then for some ρ > 0 we have

B ρ (x ) ⊂ U and max{V (x) | x ∈ B ρ (x )} = V * < h. Assume that the point x 1 is close to x such that |x 1 -x | < µ < ρ.
Without loss of generality put x 1 = 0. This particularly implies

|z 2 | < µ. (3.3)
Assume that for all x ∈ M and for all ξ ∈ R m we have

ξ T G(x)ξ ≥ c 2 |ξ| 2 . (3.4)
In addition to (1.5) assume that

V (x) ≥ c 9 , ∀x ∈ M. (3.5)
Theorem 3. If µ > 0 is small enough then the functional T attains its minimum in the set of functions

x : [s 1 , s 2 ] → U, x(•) ∈ J .
It is not hard to show that the minimum x, as long as it exists, is actually smooth. Particularly x satisfies equation (2.6).

Proof of Theorem 3. Construct a function

x + (s) = (y T + , z T + ) T (s) ∈ J as follows z + (s) = s -s 1 s 2 -s 1 z 2
and y + is determined from the Cauchy problem

(y + ) s = A(s, y + ) z 2 s 2 -s 1 , y + (s 1 ) = 0, here A(s, y) = -R -1 (x)Q(x) x=(y T ,z T + (s)) T
.

By the standard estimates we yield

|y + (s)| ≤ c 3 s -s 1 s 2 -s 1 µ, s ∈ [s 1 , s 2 ].
Here we use formula (3.3). Thus |x + (s)| < c 4 µ and for µ > 0 small enough we obtain

x + (s) ∈ B ρ (x ), ∀s ∈ [s 1 , s 2 ], x + ∈ J .
Observe that

T(x + ) ≤ c 5 µ 2 s 2 -s 1 . (3.6)
Let {u k } be a minimizing sequence for T :

T(u k ) → inf .
Then for all big enough k it follows that

T(u k ) ≤ T(x + ).
Consequently, from formulas (3.6), (3.4) and (3.2) one concludes

c 7 u k 2 C[s 1 ,s 2 ] s 2 -s 1 ≤ c 6 u k 2 ≤ µ 2 s 2 -s 1 .
This particularly implies

u k C[s 1 ,s 2 ] ≤ c 8 µ
and for µ > 0 small enough we obtain

u k (s) ∈ B ρ (x ), ∀s ∈ [s 1 , s 2 ].
The further argument is completely standard [START_REF] Taylor | Partial differential equations III[END_REF]: since the sequence

{u k } is bounded in H 1 (s 1 , s 2 ) it contains a subsequence {u k i } that is convergent in C[s 1 , s 2 ]
and weakly in H 1 (s 1 , s 2 ). Therefore the exists an element ũ ∈ H 0 such that

u k i -ũ C[s 1 ,s 2 ] → 0, s 2 s 1 ((u k i ) s -ũs ) T ϕds → 0, ∀ϕ ∈ L 2 (s 1 , s 2 ).
The element ũ is exactly the desired minimum x. This follows from the inequality:

T(u k i ) ≥ T(ũ) + 1 2 s 2 s 1 ũT s G(u k i ) h -V (u k i ) - G(ũ) h -V (ũ) ũs ds + s 2 s 1 ũT s G(u k i ) h -V (u k i ) - G(ũ) h -V (ũ) (u k i ) s -ũs ds + s 2 s 1 ũT s G(ũ) h -V (ũ) (u k i ) s -ũs ds.
The last three integrals vanish as k i → ∞. The theorem is proved.

Some Useful Facts From the Calculus of Variations

Here we collect several standard facts from the Calculus of Variations.

Let Ω ⊂ R m be an open domain with standard coordinates

x = (x 1 , . . . , x m ) T .

To proceed with formulations we split the vector x in two parts as above (2.3). Let F : Ω × R m → R be a smooth function.

We are about to state the variational problem for the functional

F x(•) = s 2 s 1 F (x(s), x s (s))ds (4.1)
with boundary conditions

z(s 1 ) = z 1 , z(s 2 ) = z 2 , y(s 1 ) = y 1 , s 1 < s 2 (4.2)
and constraints a(x, x s ) = 0. (4.3) Here a = (a 1 , . . . , a n ) T is a vector of functions that are smooth in Ω × R m .

There also must be

(y T 1 , z T 1 ) T ∈ Ω, {(y T , z T ) T | z = z 2 } ∩ Ω = ∅. Assume that det ∂a ∂y s (x, x s ) = 0, (x, x s ) ∈ Ω × R m (4.4)
and equation (4.3) can equivalently be written as

y s = Φ(y, z, z s ). Definition 1. Let a smooth function x : [s 1 , s 2 ] → Ω, x(s) = (ỹ T , zT ) T (s) be such that a(x(s), xs (s)) = 0, x(s 1 ) = x 1 = (y T 1 , z T 1 ) T , z(s 2 ) = z 2 .
We shall say that x is a stationary point of functional (4.1) with constraints (4.3) and boundary conditions (4.2) if the following holds.

For any smooth function

X : [s 1 , s 2 ] × (-ε 0 , ε 0 ) → R m , X(s, ε) = (Y T , Z T ) T (s, ε), ε 0 > 0 such that 1) X [s 1 , s 2 ] × (-ε 0 , ε 0 ) ⊂ Ω; 2) X(s, 0) = x(s), s ∈ [s 1 , s 2 ]; 3) Z(s 1 , ε) = z 1 , Z(s 2 , ε) = z 2 , Y (s 1 , ε) = y 1 , ε ∈ (-ε 0 , ε 0 ) ; 4) a(X(s, ε), X s (s, ε)) = 0, (s, ε) ∈ [s 1 , s 2 ] × (-ε 0 , ε 0 ) we have d dε ε=0 F X(•, ε) = 0.
The functions X with properties 1)-4) are referred to as variations. This theorem remains valid if the functions a, F depend on s.

For completeness of the exposition sake we prove this theorem in section 4.2.

4.1. The Homogeneous Case. In this section it is reasonable to assume the second argument of the functions a, F to be defined on a conic domain K ⊂ R m . All the formulated above results and the argument of section 4.2 remain valid under such an assumption.

Recall that by definition the domain K is a conic domain iff

x ∈ K =⇒ αx ∈ K, ∀α > 0.
Proposition 1. Assume that the function a is homogeneous in the second argument:

a(x, αx s ) = αa(a, x s ), ∀α > 0, ∀(x, x s ) ∈ Ω × K. (4.6)
Then the stationary point x preserves the "energy":

H(x, x s ) = ∂F ∂x s x s -F that is H(x(s), xs (s)) = const.
Proof of Proposition 1. Consider a function X(s, ε) = x(s + εϕ(s)) with a smooth function ϕ such that supp ϕ ⊂ [s 1 + s , s 2 -s ] and |ε|, s > 0 are small enough.

The function X satisfies all the conditions of Definition 1. Property (4.6) is needed to check condition 4) of Definition 1. Here we use integration by parts. Since ϕ is an arbitrary function the proposition is proved.

4.1.1. The Both Functions a, F are Homogeneous in x s .

Proposition 2. Let x(s) be a stationary point of F. Define a function

x(ξ) = x(f (ξ))
,

where f : [ξ 1 , ξ 2 ] → [s 1 , s 2 ] is a smooth function, f ξ (ξ) > 0, f (ξ r ) = s r , r = 1, 2.
Then x is a stationary point of F with the integral taken over [ξ 1 , ξ 2 ].

Indeed, such a reparameterization neither changes the shape of constraints (4.3) nor the shape of integral (4.1):

s 2 s 1 F (X(s, ε), X s (s, ε))ds = ξ 2 ξ 1 F ( X(ξ, ε), Xξ (ξ, ε))dξ, X(ξ, ε) = X(f (ξ), ε). Proposition 3. If F (x(s), xs (s)) > 0, s ∈ [s 1 , s 2 ]
then for any constant c > 0 we can choose a parametrization of x such that

F (x(ξ), xξ (ξ)) = c, ξ ∈ [ξ 1 , ξ 2 ].
Indeed, the desired parametrization is obtained from the equation dξ ds = 1 c F (x(s), xs (s)).

Introduce a functional

P x(•) = s 2 s 1 F (x(s), x s (s)) 2 ds
with the same constraints and boundary conditions (4.3), (4.2).

Proposition 4. Let x(s) be a stationary point of F and assume that x(s) is parametrized in accordance with Proposition 3:

F (x(s), xs (s)) = c. Then x(s) is a stationary point of P. Indeed, d dε ε=0 P X(•, ε) = 2c d dε ε=0 F X(•, ε) = 0. (4.7)
Proposition 5. Let x * (s) be a stationary point of P. Then F is the energy integral:

F (x * (s), x * s (s)) = const. Indeed, since F is homogeneous in x s it follows that ∂F ∂x s x s = F.
And the assertion follows from Proposition 1:

H = ∂F 2 ∂x s x s -F 2 = F 2 .
Proposition 6. Let x * (s) be a stationary point of P such that F (x * (s), x * s (s)) = c > 0. Then it is a stationary point of F. Indeed, it follows from (4.7). Summing up we obtain the following proposition. Proposition 7. If x(s) is a stationary point of F and F (x(s), xs (s)) > 0 then after some reparametrization it is a stationary point of P.

If x * (s) is a stationary point of P and F (x * (s), x * s (s)) = c > 0 then it is a stationary point of F. Remark 1. That is why we can not impose condition x(s 2 ) = x 2 as it is usually done for the holonomic case. The value Y (s 2 , ε) has already been uniquely defined by other boundary conditions and the constraints.

In other words if we add the condition Y (s 2 , ε) = y 2 to the conditions 1)-4) of Definition 1 then the set of variations {X(s, ε)} may turn up to be insufficiently large to prove theorem 4.

For example, consider a plane R 2 = {x = (y, z) T }. There is a unique smooth path from x 1 = (0, 0) T to x 2 = (0, 1) T that satisfies the equation y s = 0. (Much more complicated example by C. Caratheodory see in [START_REF] Arnold | Neishtadt: Dynamical Systems III[END_REF].) Cauchy problem (4.9) has the suitable solution at least for |ε| and s 2 -s 1 small. Observe also that Y ε (s 1 , ε) = 0. This is a system of linear ordinary differential equations for λ. Due to assumption (4.4) this system can be presented in the normal form that is λ s = Λ(s, λ). Since we know λ(s 2 ), by the existence and uniqueness theorem we obtain λ(s) as a solution to the IVP for (4.14).

Equation (4.13) takes the form

s 2 s 1 [F * ] z δzds = 0.
Since δz is an arbitrary function we get [F * ] z = 0. Together with (4.14) this proves the theorem.

2 s 1 H

 21 Furthermore we haveX = x(s) + εϕ(s)x s (s) + O(ε 2 ), X s = xs (s) + ε ϕ s (s)x s (s) + ϕ(s)x ss (s) + O(ε 2 ) xs ) + ϕ s (s) ∂F (x,xs ) ∂x s xs (s) ds = s (x(s), xs (s))ϕ s (s)ds = 0.

4. 2 .

 2 Proof of Theorem 4. Introduce a notation [F ] y = -F ] x = ([F ] y , [F ] z ). Let us put Z(s, ε) = z(s) + εδz(s), supp δz ⊂ [s 1 , s 2 ]. (4.8)Then the function Y is uniquely determined from the following Cauchy problemY s (s, ε) = Φ(Y (s, ε), Z(s, ε), Z s (s, ε)), Y (s 1 , ε) = y 1 .(4.9)

(4. 10 ) 1 [ 2 s 1 λdε ε=0 A = s 2 s 1 [ 2 s 1 [F

 101212121 Using the standard integration by parts technique and from formulas(F ] z δz + [F ] y Y ε ds + ∂F ∂y s (x(s 2 ), xs (s 2 ))Y ε (s 2 , 0) = 0. (4.11)The function λ(s) is still undefined but due to condition (4.4) the value λ(s 2 ) is determined uniquely from (4.5).From condition 4) of definition 1 it follows thatA(ε) = s (s)a(X(s, ε), X s (s, ε))ds = 0.By the same argument as above we haved λa] z δz + [λa] y Y ε ds + λ(s 2 ) ∂a ∂y s (x(s 2 ), xs (s 2 ))Y ε (s 2 , 0) = 0. (4.12)Summing formulas (4.12) and (4.11) we yields * ] z δz + [F * ] y Y ε ds = 0. (4.13) To construct the function λ consider an equation [F * ] y = 0. (4.14)

  Theorem of Existence of Minimum. Let U ⊂ R m be a coordinate patch in M with coordinates (2.3). Let | • | stand for the standard Euclidean norm in R m and let B r (x) ⊂ R m stand for the open ball of radius r and with center in x.

1 , s 2 ) equation (2.2) holds. 3.2.

  (s, x, x s ) = F (x, x s ) + λ(s)a(x, x s ),

	Theorem 4 ([2]). If the function x is a stationary point of functional
	(4.1) with constraints (4.3) and boundary conditions (4.2) then there
	is a smooth function λ(s) = (λ 1 , . . . , λ n )(s) such that x satisfies the
	equations
	d ds = 0, F and ∂F * ∂x s -∂F * ∂x ∂F ∂y s (x(s 2 ), xs (s 2 )) + λ(s 2 )	∂a ∂y s	(x(s 2 ), xs (s 2 )) = 0.	(4.5)

*

The research was funded by a grant from the Russian Science Foundation (Project No. 19-71-30012).


Oleg Zubelevich
email: ozubel@yandex.ru

A Note on the Variational Problem of Lagrange

It is shown that the Lagrange problem with differential constraints on an m-dimensional manifold is reduced to a Hamiltonian system on a 2m-dimensional symplectic manifold.

Introduction. The Statement of the Problem

Let M be a smooth m-dimensional manifold with local coordinates x = (x 1 , . . . , x m ) T . We assume all the objects below to be smooth enough for the formulas to make sense.

Let L : T M → R, L = L(x, v) be a Lagrangian function. Assume that the mapping

(x, v) → (x, y), y = ∂L ∂v (x, v) (1.1)
is a diffeomorphism of T M to T * M and equation (1.1) implies v = f (x, y), f = (f 1 , . . . , f m ) T .

Here (x, y) are the local coordinates in T * M .

2010 Mathematics Subject Classification. 70H05, 70H03, 34C40, 53D05, 34B15. Key words and phrases. Lagrange problem, variational calculus, symplectic geometry, Hamiltonian systems .

The research was funded by a grant from the Russian Science Foundation (Project No. 19-71-30012).

Assume also that

A(x, v) = ∂ 2 L ∂v 2 (x, v) is a matrix of positively definite quadratic form; ∂f ∂y = A -1 . (1.2) Introduce 1-forms b i = b i j (x)dx j , i = 1, . . . , n < m; and B(x) =     b 1 1 (x) b 1 2 (x) • • • b 1 m (x) b 2 1 (x) b 2 2 (x) • • • b 2 m (x) . . . . . . . . . . . . b n 1 (x) b n 2 (x) • • • b n m (x)     , rang B(x) = n, x ∈ M.
These forms define an m -n-dimensional differential system

T x = n i=1 ker b i (x) ⊂ T x M in M [3].
Let Σ ⊂ M be a smooth n-dimensional submanifold such that

x ∈ Σ =⇒ T x M = T x Σ ⊕ T x .
Let S stand for a set of functions x : [t 1 , t 2 ] → M such that

x(t 1 ) = x 1 , x(t 2 ) ∈ Σ, ẋ(t) ∈ T x(t) , t ∈ [t 1 , t 2 ].
Here x 1 ∈ M is a fixed point.

In coordinate notation the last inclusion is read as follows

B x(t) ẋ(t) = 0. (1.3)
The scheme described below works if we replace constraints (1.3) with ones nonlinear in ẋ. We use linear constraints just to avoid an irrelevant complication of the technique.

From the calculus of variations [START_REF] Akhiezer | The Calculus of Variations[END_REF] we know that if a function x(•) ∈ S is a critical point of the functional

F : S → R, F x(•) = t 2 t 1 L x(t), ẋ(t) dt then there is a vector λ(t) = (λ 1 , . . . , λ n )(t) such that the function x(t) satisfies the Lagrange equations d dt ∂L * ∂v x, ẋ, λ(t) - ∂L * ∂x x, ẋ, λ(t) = 0, (1.4) with L * (x, v, λ) = L(x, v) + λB(x)v, L * : T M × R n → R; and ∂L * ∂v x(t 2 ), ẋ(t 2 ), λ(t 2 ) ξ = 0, ξ ∈ T x(t 2 ) Σ. (1.5)
The main object of our study is system (1.3), (1.4). Loosely speaking our aim is to show that this system is described in terms of a Hamiltonian system in some 2m-dimensional symplectic manifold.

As far as the author knows this simple but important fact has not been noticed anywhere.

Note also that the variational problem described above provoked a lot of mistakes in classical mechanics. Many researches thought that such a variational problem is equivalent to the Lagrange-d'Alembert equations for the mechanical system with the Lagrangian L and ideal constraints (1.3). In general it is not so. This topic together with nice historical details is discussed in [START_REF] Bloch | Nonholonomic Mechanics and Control (Interdisciplinary Applied Mathematics)[END_REF].

The Main Construction

Let H * : Q → R, Q = T * M × R n
stand for the standard Legendre transformation of L * :

p = (p 1 , . . . , p m ) = ∂L * ∂v (x, v, λ), v = f (x, p, λ), where f (x, p, λ) := f x, p -λB(x) ; H * (x, p, λ) = p f (x, p, λ) -L x, f (x, p, λ) -λB(x) f (x, p, λ).
The following well-known fact holds true. Let the function λ = λ(t) be fixed then the function x(t) is a solution to equations (1.4) iff the functions

x(t), p(t) = ∂L * ∂ ẋ x(t), ẋ(t), λ(t)
are the solution to the Hamilton equations ẋ = ∂H * ∂p x, p, λ(t) , ṗ = -∂H * ∂x x, p, λ(t) .

(2.1)

Remark 1. The condition (1.5) takes the form

p(t 2 )ξ = 0, ξ ∈ T x(t 2 ) Σ.
Equation (1.3) takes the form

B(x) f (x, p, λ) = 0. (2.2)
To determine the function λ(t) one must differentiate this equation in t and substitute ẋ, ṗ from (2.1). This procedure gives an equation of the form

R(x, p, λ) λT = ψ(x, p, λ), R(x, p, λ) = B(x) ∂f ∂y x, p -λB(x) B T (x).
Recall that A is a matrix of positive definite quadratic form. So that by formula (1.2) we have det R(x, p, λ) = 0;

(2.3) and the last equation takes the form λ = Λ(x, p, λ).

(2.4) System (2.1), (2.4) is an ODE system on the manifold Q with the local coordinates (x, p, λ).

Observe that the way we construct system (2.1), (2.4) implies that the vector function

C(x, p, λ) = B(x) f (x, p, λ)
is a first integral of this system. This in particular means that a solution (x, p, λ)(t) to system (2.1), (2.4) satisfies (2.2) identically provided the solution's initial conditions satisfy (2.2).

Let N denote the zero level manifold of C(x):

N = {(x, p, λ) ∈ Q | C(x, p, λ) = 0}.
Proposition 1. The set N is a smooth 2m-dimensional submanifold of Q; this submanifold is locally presented as a graph of a function: λ = ϕ(x, p).

(2.5) Indeed, this follows from the Implicit Function Theorem and formula (2.3) : ∂C ∂λ = R.

Remark 2. Formula (2.5) implies that (x, p) are the local coordinates in N .

The manifold T * M has the standard symplectic structure:

Ω = dx i ∧ dp i .
This 2-form is trivially extended to the cross product Q. Therefore we proceed with considering of Ω in Q.

Proposition 2. The manifold N is a symplectic manifold with a symplectic form given by the following restriction

ω = Ω | N .
The coordinates (x, p) are symplectic coordinates in N .

Indeed, this is a direct consequence of proposition 1 and remark 2. Let w(z) ∈ T z Q, z = (x, p, λ) denote the vector field of system (2.1), (2.4). Since C is a first integral we have z ∈ N ⇒ w(z) ∈ T z N . Theorem 1. The vector field u = w | N is a Hamiltonian vector field in N . The Hamiltonian of u is

H = H * | N : N → R.
In the local coordinates (x, p) this Hamiltonian is given by the formula H(x, p) = H * x, p, ϕ(x, p) .

Proof of Theorem 1. Let i w stand for the interior product. The following formula is the main: The theorem is proved.

Example

Let the Lagrangian L be given by the formula

L(x, ẋ) = 1 2 ẋT G(x) ẋ,
where G(x) is a Riemann metric in M . After some calculation one yields H(x, p) = 1 2 pP(x)G -1 (x)p T ,

  i w Ω = dH * -∂H * ∂λ j dλ j . (2.6) This formula is checked by direct calculation. The formulas ∂H * ∂λ = -B(x) f (x, p, λ), ∂H * ∂λ N = 0 are evident. Restricting the both sides of formula (2.6) to N we obtain i u ω = dH.

where

is the orthogonal projection: