Jan Soko Lowski
email: jan.sokolowski@iecn.u-nancy.fr

Jan Stebel
email: stebel@math.cas.cz

Shape sensitivity analysis of incompressible non-Newtonian fluids

Keywords: shape optimization, shape gradient, incompressible fluid, non-Newtonian fluid, Navier-Stokes equations

We study the shape differentiability of a cost function for the steady flow of an incompressible viscous fluid of power-law type. The fluid is confined to a bounded planar domain surrounding an obstacle. For smooth perturbations of the shape of the obstacle we express the shape gradient of the cost function which can be subsequently used to improve the initial design.

Introduction

Shape optimization for nonlinear partial differential equations is a growing field in the contemporary optimum design of structures. In this field systems of the solid and fluid mechanics as well as e.g., the coupled models of fluid-structure interaction are included for real life problems. The main difficulty associated with the mathematical analysis of nonlinear state equations is the lack of existence of global strong solutions for mathematical models in three spatial dimensions.

In numerical methods of shape optimization the common approach is the discretization of continuous shape gradient. Therefore, the proper derivation and analysis of the regularity properties of the shape gradient is crucial for numerical solution of the shape optimization problem. The shape sensitivity analysis requires, in particular, the proof of the Lipschitz continuity of solutions the the state equations with respect to the boundary variations. This property of solutions can be obtained e.g. by analysis of the state equation transported to the fixed reference domain which is explained in the case of linear elliptic boundary value problems in monograph [START_REF] Soko | Introduction to shape optimization. Shape sensitivity analysis[END_REF]. For the nonlinear problems the Lipschitz continuity is not obvious and it requires the additional regularity of solutions to the state equation. In addition, for the applications of levelset method of shape optimization it is required that the obtained shape gradient of the cost functional is given by a function while the general theory gives only the existence of a distribution. In conlusion, it seems that the shape sensitivity analysis in the case of a nonlinear state equation is the main step towards the numerical solution of the shape optimization problems.

In gas dynamics described by the compressible Navier-Stokes there is the existence of weak global solutions. However, the shape sensitivity analysis can be performed only for specific local solutions. The state of art in shape optimization for compressible Navier-Stokes equations is presented in the forthcoming monograph [START_REF] Plotnikov | Compressible Navier-Stokes equations. Theory and shape optimization[END_REF], see also [START_REF] Plotnikov | Shape derivative of drag functional[END_REF]. For incompressible Navier-Stokes equations, the sensitivity analysis of shape functionals is performed e.g. in [START_REF] Consiglieri | New approach to the incompressible Maxwell-Boussinesq approximation: existence, uniquness and shape sensitivity[END_REF] and [START_REF] Moubachir | Moving shape analysis and control[END_REF]. In this paper we are concerned with the non-Newtonian model where the stress is a (nonlinear) function of the velocity gradient. Optimal control problem for this model was studied in [START_REF] Slawig | Distributed control for a class of non-newtonian fluids[END_REF][START_REF] Wachsmuth | Optimal control of planar flow of incompressible non-Newtonian fluids[END_REF]. Numerical shape optimization was done in [START_REF] Abraham | Shape optimization in steady blood flow: A numerical study of non-newtonian effects[END_REF], see also [START_REF] Haslinger | Shape optimization for navier-stokes equations with algebraic turbulence model: Numerical analysis and computation[END_REF]. We present new results on the existence of the shape gradient.

We consider the steady flow of an incompressible fluid in a bounded domain

Ω := B \ S in R 2
, where B is a container and S is an obstacle. Motion of the fluid is described by the system of equations

div (v ⊗ v) -div S(Dv) + ∇p + Cv = f in Ω, div v = 0 in Ω, (P (Ω)) v = g on ∂Ω.
Here v, p, C, f stands for the velocity, the pressure, the constant skew-symmetric Coriolis tensor and the body force, respectively. The traceless part S of the Cauchy stress can depend on the symmetric part Dv of the velocity gradient in the following way:

S(Dv) = ν(|Dv| 2 )Dv, (1)
where ν, |Dv| 2 is the viscosity and the shear rate, respectively. In particular, we assume that ν has a polynomial growth (see Section 2.1 below), which includes e.g. the Carreau and the power-law model.

In the model the term of Coriolis type is present. This term appears e.g. when the change of variables is performed in order to take into account the flight scenario of the obstacle in the fluid.

The aim of this paper is to investigate differentiability of a shape functional depending on the solution to (P (Ω)) with respect to the variations of the shape of the obstacle. We consider a model problem with the drag functional

J(Ω) := ∂S (S(Dv) -pI)n • d, (2)
with a given constant unit vector d. Instead of J one could take other type of functional, since our method does not rely on its specific form.

Our main interest is the rigorous analysis of the shape differentiability for (P (Ω)) and [START_REF] Consiglieri | New approach to the incompressible Maxwell-Boussinesq approximation: existence, uniquness and shape sensitivity[END_REF]. We follow the general framework developed in [START_REF] Soko | Introduction to shape optimization. Shape sensitivity analysis[END_REF] using the speed method and the notion of the material derivative. Let us point out that due to (1) the state problem is nonlinear in its nature. We refer the reader to [START_REF] Tröltzsch | Optimal control of partial differential equations[END_REF] for an introduction to optimization problems for nonlinear partial differential equations.

Shape derivatives

We start by the description of the framework for the shape sensitivity analysis. For this reason, we introduce a vector field T ∈ C 2 (R 2 , R 2 ) vanishing in the vicinity of ∂B and define the mapping

y(x) = x + εT(x).
For small ε > 0 the mapping x → y(x) takes diffeomorphically the region Ω onto Ω ε = B \ S ε where S ε = y(S). We consider the counterpart of problem (P (Ω)) in Ω ε , with the data f |Ωε and g |Ωε . The new problem will be denoted by (P (Ω ε )) and its solution by (v ε , pε ).

For the nonlinear system (P (Ω)) we introduce the shape derivatives of solutions. To this end we need the linearized system of the form:

Find the couple (u, π) such that

div [u ⊗ v + v ⊗ u -S ′ (Dv)Du] + ∇π + Cu = F in Ω, div u = 0 in Ω, (P lin (Ω)) u = h on ∂Ω,
where F and h are given elements.

The shape derivative v ′ and the material derivative v of solutions are formally introduced by

v ′ := lim ε→0 vε -v ε , v := lim ε→0 vε • y -v ε ,
where vε • y(x) := vε (y(x)), and are related to each other as follows:

v = v ′ + (∇v)T.
The standard calculus for differentiating with respect to shape yields that v ′ is the solution of (P lin (Ω)) with the data F = 0 and h = -∂v/∂n(T • n). Using [START_REF] Plotnikov | Shape derivative of drag functional[END_REF] as the definition of J we obtain the expression for the shape gradient:

dJ(Ω; T) := lim ε→0 J(Ω ε ) -J(Ω) ε = Ω [(Cv ′ ) • ξ + (S ′ (Dv)Dv ′ -v ′ ⊗ v -v ⊗ v ′ ) : ∇ξ] - ∂S (f • d)T • n. (3)
In the above formula, the part containing v ′ depends implicitly on the direction T. This is not convenient for practical use, hence we introduce the adjoint problem for further simplification of (3): Find the couple (w, s) such that

-2(Dw)v -div S ′ (Dv) ⊤ Dw + ∇s -Cw = 0 in Ω, div w = 0 in Ω, w = d on ∂Ω. (P adj (Ω))
Consequently, the expression for dJ reduces to

dJ(Ω; T) = - ∂S S ′ (Dv) ⊤ Dw -sI : ∂v ∂n ⊗ n + f • d T • n. (4)
In order to prove the result given by ( 3) and ( 4) we need the material derivatives. In particular, it is sufficient to show that the linear mapping

T → dJ(Ω; T)
is continuous in an appropriate topology, see the structure Theorem in the book [START_REF] Soko | Introduction to shape optimization. Shape sensitivity analysis[END_REF] for details.

Preliminaries

We impose the structural assumptions on the data, state the known results on well-posedness of (P (Ω)) and introduce the elementary notation for shape sensitivity analysis.

Structural assumptions

We require that S has a potential

Φ : [0, ∞) → [0, ∞), i.e. S ij (D) = ∂Φ(|D| 2 )/∂D ij .
Further we assume that Φ is a C 3 function with Φ(0) = 0 and that there exist constants C 1 , C 2 , C 3 > 0 and r ≥ 2 such that

C 1 (1 + |A| r-2 )|B| 2 ≤ S ′ (A) :: (B ⊗ B) ≤ C 2 (1 + |A| r-2 )|B| 2 , ( 5a
)
|S ′′ (A)| ≤ C 3 (1 + |A| r-3 ) (5b)
for any 0 = A, B ∈ R 2×2 sym . Here the symbol :: stands for the usual scalar product in R 2 4 . The above inequalities imply the monotone structure of S, see e.g. [START_REF] Málek | Mathematical issues concerning the Navier-Stokes equations and some of its generalizations[END_REF].

Weak formulation

For the definition of the weak solution we will use the space

W 1,r 0,div (Ω) := {φ ∈ W 1,r 0 (Ω); div φ = 0}.
Let f ∈ (W 1,2 0,div (Ω)) * and g ∈ W 1,r (Ω) with div g = 0. Then a function v ∈ g + W 1,r 0,div (Ω) is said to be a weak solution to the problem (P (Ω)) if

Ω S(Dv) : Dφ -v ⊗ v : ∇φ + Cv • φ = Ω f • φ (6)
for every φ ∈ W 1,r 0,div (Ω). Note that the pressure is eliminated since test functions are divergence free.

The following result was shown in [START_REF] Kaplický | On global existence of smooth twodimensional steady flows for a class of non-newtonian fluids under various boundary conditions[END_REF].

Theorem 1 (Kaplický et al. [START_REF] Kaplický | On global existence of smooth twodimensional steady flows for a class of non-newtonian fluids under various boundary conditions[END_REF]). Let Ω ∈ C 2 , f ∈ L 2+ǫ0 (Ω), ǫ 0 > 0 and (5a)-(5b) hold with r > 3 2 . Then there exists a constant δ > 0 such that for every g satisfying g 3,q ≤ δ, (q > 2), problem (P (Ω)) has a weak solution satisfying v ∈ W 2,2+ǫ (Ω), p ∈ W 1,2+ǫ (Ω), ǫ > 0.

Note that the above result applies only to the unperturbed domain, i.e. ε = 0. Assuming smallness of f 2,B and g 3,q,B , one can prove that (P (Ω)), (P (Ω ε ) has a unique weak solution satisfying

v ≤ C E ( f 2,B , g 3,q,B ) and vε ≤ C E ( f 2,B , g 3,q,B ),
respectively, where C E is independent of ε. At this point we summarize the main hypotheses.

Assumption 1

In what follows, Ω ∈ C 2 is a bounded planar domain of the form

Ω = B \ S, f ∈ L 2+ǫ0 (B), ǫ 0 > 0, g ∈ W 3,q (B) (q > 2)
is supported in the vicinity of ∂B, (5a)-(5b) hold with r ∈ [2, 4) and f 2,B , g 3,q,B are small enough.

Let us point out that equation ( 2) which defines J is not suitable for weak solutions in general, since the energy inequality does not provide enough information about the trace of p and Dv. We therefore introduce an alternative definition that requires less regularity. Let us fix an arbitrary divergence free function ξ ∈ C ∞ c (B, R 2 ) such that ξ = d in a vicinity of S. Then, integrating (2) by parts and using (P (Ω)) yields:

J(Ω) = Ω [(Cv -f ) • ξ + (S(Dv) -v ⊗ v) : ∇ξ] . (7)
Note that this identity is finite for any v ∈ W 1,2 (Ω).

Deformation of the shape

Let us introduce the following notation: We will denote by D T the Jacobian matrix whose components are

(D T) ij = (∇T) ji = ∂ i T j . Further, N(x) := g(x)M -1 (x), M(x) := I + ε D T(x), g(x) := det M(x).
One can easily check that the matrix N and the determinant g admit the expansions:

g = 1 + εdiv T + O(ε 2 ), N = I + εN ′ + O(ε 2 ), N ′ = (div T)I -D T, (8)
where the symbol O(ε 2 ) denotes a function whose norm in C 1 (Ω) is bounded by Cε 2 , see [START_REF] Soko | Introduction to shape optimization. Shape sensitivity analysis[END_REF].

Theorem 3. Let Assumption 1 be satisfied and f ∈ W 1,2 (Ω). Then the shape gradient of J reads dJ(Ω, T) = J v (ṽ) + J e (T),

where the dynamical part J v and the geometrical part J e is given by

J v (ṽ) = Ω [(Cṽ) • ξ + (S ′ (Dv)Dṽ -ṽ ⊗ v -v ⊗ ṽ) : ∇ξ] , J e (T) = Ω (I tr N ′ -N ′ ) Cv -CN ′⊤ v -(I tr N ′ -N ′ ) f -(∇f )T • ξ + v ⊗ N ′⊤ v + S ′ (Dv) (N ′ ∇v -∇(N ′⊤ v)) sym -(tr N ′ )Dv + N ′⊤ S(Dv) : ∇ξ + v ⊗ v -S(Dv) : ∇(N ′⊤ ξ) ,
respectively. In particular, as ṽ depends continuously on T, the mapping

T → dJ(Ω, T) is a bounded linear functional on C 2 (R 2 , R 2 ).
Based on the previous result we can deduce that the shape gradient has the form of a distribution supported on the boundary of the obstacle. Since this representation is unique, the formal results derived in Section 1.1 are justified provided that the shape derivatives and adjoints exist and are sufficiently regular.

Corollary 1. Let Assumption 1 be satisfied. Then (i) the shape derivative v ′ exists and is the unique weak solution to (P lin (Ω))

with F = 0, h = - ∂v ∂n (T • n);
(ii) the adjoint problem (P adj (Ω)) has a unique weak solution that satisfies: w ∈ W 2,2 (Ω) and s ∈ W 1,2 (Ω).

If in addition f ∈ W 1,2 (Ω), then (iii) the shape gradient of J satisfies (3);

(iv) the representation (4) is satisfied in the following sense:

dJ(Ω; T) = - ∂S S ′ (Dv) ⊤ Dw -sI : ∂v ∂n ⊗ n + f • d T • n. (12)
In the remaining part we show the main steps of the proof of Theorem 3. Details can be found in [START_REF] Soko | Shape sensitivity analysis of time-dependent flows of incompressible non-Newtonian fluids[END_REF], where the time-dependent problem is treated.

Formulation in the fixed domain

In this section we transform the problem (P (Ω ε )) to the fixed domain Ω. Let us introduce the following notation:

v ε (x) := N ⊤ (x)v ε (y(x)), x ∈ Ω.
Note that the definition of v ε implies that div v ε = 0. The new function v ε ∈ g + W 1,r 0,div (Ω) satisfies the equality

Ω gS(D ε v ε ) : D ε φ -v ε ⊗ v ε : ∇φ + Cv ε • φ = Ω f • φ + A 1 ε , φ W 1,2 0,div (Ω) for all φ ∈ W 1,r 0,div (Ω), ( 13
)
where the term A 1 ε on the right hand side is defined for φ ∈ W 1,2 0,div (Ω) by

A 1 ε , φ W 1,2 0,div (Ω) = Ω v ε ⊗ N -⊤ v ε : ∇(N -⊤ φ) -v ε ⊗ v ε : ∇φ + (C -gN -1 CN -⊤ )v ε • φ + (gN -1 f • y -f ) • φ . (14)
Here

D ε v ε := g -1 (N∇(N -⊤ v ε )) sym .
Applying change of coordinates we further get:

J(Ω ε ) = Ω g N -1 CN -⊤ v ε -N -1 f • y • ξ + N ⊤ S(D ε v ε ) -v ε ⊗ (N -⊤ v ε ) : ∇(N -⊤ ξ) . ( 15
)
Now after all quantities and equations have been transformed to the fixed domain Ω, we can analyze the limit ε → 0. Lemma 1. The sequence {v ε } ε>0 is bounded in W 1,r 0,div (Ω) and satisfies:

v ε ⇀ v weakly in W 1,r 0,div (Ω), N ⊤ S(D ε v ε ) ⇀ S(Dv) weakly in L r ′ (Ω, R 2×2 ), A 1 ε ⇀ 0 weakly in W 1,r 0,div (Ω) * .
In particular, v is the unique weak solution to (P (Ω)).

Existence of material derivative

Our next task is to identify ṽ as the limit of the sequence {u ε }, where

u ε := v ε -v ε .
First we write down the system for the differences u ε . Subtracting ( 13) and ( 6) we find that u ε ∈ W 1,r 0,div (Ω) satisfies the equality

Ω 1 ε g(S(D ε v ε ) -S(D ε v)) : D ε φ + Cu ε • φ -(v ε ⊗ u ε + u ε ⊗ v) : ∇φ = 1 ε A ε , φ W 1,2 0,div (Ω) (16)
for all φ ∈ W 1,r 0,div (Ω). The term A ε ∈ W 1,2 0,div (Ω) * on the right hand side is defined as follows:

A ε := A 1 ε + A 2 ε , A 1
ε is given by ( 14), Next we state the properties of the sequence {u ε } ε>0 .

A 2 ε , φ W 1,
Lemma 2. The sequence {u ε } ε>0 is bounded in W 1,2 0,div (Ω). Further it holds:

A ε ε ⇀ A ′ 0
weakly in W 1,2 0,div (Ω) * , u ε ⇀ ṽ weakly in W 1,2 0,div (Ω), 1 ε (g(S(D ε v ε ) -S(D ε v)), D ε φ) → (S ′ (Dv)Dṽ, Dφ) for all φ ∈ W 1, 2r 4-r (Ω),

where A ′ 0 is defined in (10a) and ṽ is the solution of (P lin (Ω)) with F := A ′ 0 and h = 0.

This completes the proof of Theorem 2.

Shape gradient of J

To prove Theorem 3, we decompose the fraction

J(Ω ε ) -J(Ω) ε = J ε 1 + J ε 2
in a suitable way. Using Lemma 1 and Lemma 2 and the properties of g and N ′ , it is then possible to show that J ε 1 → J v (ṽ) and J ε 2 → J e (T). The continuity of the map T → dJ(Ω; T) follows from the estimate [START_REF] Soko | Introduction to shape optimization. Shape sensitivity analysis[END_REF].

  2 0,div (Ω) := Ω N ⊤ S(D ε v) : ∇(N -⊤ φ) -S(Dv) : Dφ .

Acknowledgement

The work of J. Stebel was supported by the Czech Science Foundation (GA ČR) grant No. 201/09/0917 and by the ESF grant Optimization with PDE Constraints. The research of J. Soko lowski is performed in the framework of the ANR project GAOS.

The value of the shape functional for Ω ε is given by

where ξ ε := (N -⊤ ξ) • y -1 . Using the properties of the Piola transform one can check that div ξ ε = 0. If vε and pε were sufficiently smooth, it would hold that

Nevertheless, as opposed to (P (Ω)), we do not require any additional regularity of the solution to the perturbed problem (P (Ω ε )) and hence the expression in [START_REF] Slawig | Distributed control for a class of non-newtonian fluids[END_REF] need not be well defined. We introduce the auxiliary function ṽ:

which is related to the material derivative v by the identity

For the justification of the results of the paper we will use ṽ since, unlike the material derivative, it preserves the divergence free condition.

Main results

The first result is the existence of ṽ and hence also of the material derivative.

Theorem 2. Let Assumption 1 be satisfied. Then the function ṽ exists and is the unique weak solution of (P lin (Ω)) with the data

The following estimate holds:

Next we establish the existence of the shape gradient of J.


Tadeusz Antczak

The Exact l 1 Penalty Function Method for Constrained Nonsmooth Invex Optimization Problems

Keywords: exact l1 penalty function method, absolute value penalty function, penalized optimization problem, locally Lipschitz invex function, Generalized Karush-Kuhn-Tucker optimality conditions

The exactness of the penalization for the exact l1 penalty function method used for solving nonsmooth constrained optimization problems with both inequality and equality constraints are presented. Thus, the equivalence between the sets of optimal solutions in the nonsmooth constrained optimization problem and its associated penalized optimization problem with the exact l1 penalty function is established under locally Lipschitz invexity assumptions imposed on the involved functions.

Introduction

Considerable attention has been given in recent years to devising methods for solving nonlinear programming problems via unconstrained minimization techniques. One class of methods which has emerged as very promising is the class of exact penalty function methods. Methods using exact penalty function transform a constrained extremum problem into a single unconstrained optimization problem. The constraints are placed into the objective function via a penalty parameter c in a way that penalizes any violation of the constraints.

One important property that distinguishes exact penalty functions is the exactness of the penalization. The concept of exact penalization is sometimes ambiguous, or at least varies from author to author. One of the de…nitions of the exactness of the penalization is the following: there is an appropriate penalty parameter choice such that a single unconstrained minimization of the penalty function yields a solution of the constrained optimization problem.

Nondi¤erentiable exact penalty functions were introduced for the …rst time by Eremin [START_REF] Eremin | The penalty method in convex programming[END_REF] and Zangwill [START_REF] Zangwill | Nonlinear programming via penalty functions[END_REF]. In almost all of the introduced penalized approaches the notion of convexity plays a dominant role. In 1970, Luenberger [START_REF] Luenberger | Control problem with kinds[END_REF] showed that, under convex assumptions, there is a lower bound for a penalty parameter c, equal to the largest Lagrange multiplier in absolute value, associated to the nonlinear optimization problem. Later, Charalambous [START_REF] Charalambous | A lower bound for the controlling parameters of the exact penalty functions[END_REF] generalized the result of Luenberger for the absolute value penalty function, assuming the second-order su¢cient conditions. Under the assumptions that the minimization problem is solvable and that it satis…es the relaxed Slater constraint qual-i…cation, Mangasarian [START_REF] Mangasarian | Su¢ciency of exact penalty minimization[END_REF] characterized solutions of the convex optimization problem in terms of minimizers of the exact penalty function for a single value of the penalty parameter exceeding some treshold. Bazaraa et al. [START_REF] Bazaraa | Nonlinear programming: theory and algorithms[END_REF] also used the exact l 1 penalty function method to solve nonlinear convex optimization problems with both inequality and equality constraints. They assumed that the objective function and the inequality constraints are convex and the equality constraints are a¢ne functions to prove that a Karush-Kuhn-Tucker point in the original optimization problem is a minimizer of the exact l 1 penalty function in the associated penalized optimization problem with su¢ciently large value of a penalty parameter. In the mentioned above works, the lower bound of the penalty parameter above which, for all penalty parameters, any optimal solution of the original nonlinear optimization problem is also a minimizer of the penalized problem has been given for di¤erentiable optimization problems involving convex functions. However, from the practical point of view, the converse result is also important.

In recent years, some numerous generalizations of convex functions have been derived which proved to be useful for extending optimality conditions and some classical duality results, previously restricted to convex programs, to larger classes of nonconvex optimization problems. One of them is the invexity notion introduced by Hanson [START_REF] Hanson | On su¢ciency of the Kuhn-Tucker conditions[END_REF] for di¤erentiable scalar functions and later generalized from di¤erent points of view, also in the case of nondi¤erentiable functions (see [START_REF] Antczak | Lipschitz r-invex functions and nonsmooth programming[END_REF], [START_REF] Craven | Invex functions and constrained local minima[END_REF], [START_REF] Egudo | On su¢ciency of Kuhn-Tucker conditions in nonsmooth multiobjective programming[END_REF], [START_REF] Kaul | Generalized nonsmooth invexity[END_REF], [START_REF] Kim | On duality theorems for nonsmooth Lipschitz optimization problems[END_REF], [START_REF] Reiland | Nonsmooth invexity[END_REF], [START_REF] Zhao | On su¢ciency of the Kuhn-Tucker conditions in nondi¤erentiable programming[END_REF], and others). Now, we show that there is the equivalence between the set of optimal solutions in a nondi¤erentiable nonconvex optimization problem and the set of minimizers in its associated exact penalized problem with the absolute value penalty function. It turns out that this property is not true only for (di¤erentiable) convex optimization problems, but it still holds for nonlinear optimization problems involving locally Lipschitz invex functions with respect to the same function (with the exception of those equality constraint functions for which the associated Lagrange multipliers are negative -these functions should be assumed to be incave with respect to the same function ). The result established here shows that there does exist a lower bound for a penalty parameter c, equal to the largest Lagrange multiplier in absolute value, associated to a Karush-Kuhn-Tucker point in the original nonlinear optimization problem, above which this equivalence holds. Further, in the case when at least one of the functions constituting the nondi¤erentiable constrained optimization problem is not locally Lipschitz invex and in the case when the objective function is coercive but not invex, then the equivalence in the sense discussed here might not hold between these optimization problems.

Preliminaries and Problem Formulation

Throughout this section, X is a nonempty subset of R n . A real-valued function f : X ! R is said to be locally Lipschitz on X if, for any x 2 X, there exist a neighborhood U of x and a positive constant K x > 0 such that, for every y; z 2 U , it holds jf (y) f (z)j 5 K x ky zk. The Clarke generalized directional derivative [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF] of a locally Lipschitz function f :

X ! R at x 2 X in the direction v 2 R n , denoted f 0 (x; v), is given by f 0 (x; v) = lim sup y!x #0 f (y+ v) f (y) .
De…nition 1 The Clarke generalized subgradient [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF]

of f at x 2 X, denoted @f (x), is de…ned by @f (x) = n 2 R n : f 0 (x; v) T v for all v 2 R n o .
The following de…nition is a generalization of the de…nition of a class of di¤erentiable convex functions to the case of a class of locally Lipschitz invex functions (see [START_REF] Hanson | On su¢ciency of the Kuhn-Tucker conditions[END_REF]).

De…nition 2 [START_REF] Hanson | On su¢ciency of the Kuhn-Tucker conditions[END_REF] Let a function f : X ! R be a locally Lipschitz function on X and u 2 X. If there exists a vector-valued function : X X ! R n such that, for each x 2 X, the inequality f (x) f (u)

T (x; u) holds for any 2 @f (u), then f is said to be a locally Lipschitz invex function at u on X with respect to . If the inequality above is satis…ed at any point u, then f is said to be a locally Lipschitz invex function on X with respect to .

In order to de…ne an analogous class of (strictly) Lipschitz incave functions with respect to , the direction of the inequality in the de…nition of invex functions should be changed to the opposite one.

De…nition 3 [15] A continuous function

f : R n ! R is said to be coercive if lim kxk!1 f (x) = 1:
Consider the following constrained optimization problem: minimize f (x) subject to g i (x) 0, i 2 I = f1; :::mg ; h j (x) = 0, j 2 J = f1; :::; sg x 2 X;

(P)
where f : X ! R and g i : X ! R, i 2 I, h j : X ! R, j 2 J, are locally Lipschitz functions on a nonempty set X R n . Let D := fx 2 X : g i (x) 0, i 2 I, h j (x) = 0, j 2 Jg be the set of all feasible solutions of problem (P). Further, we denote a set of active inequality constraints at point x 2 X by I (x) = fi 2 I : g i (x) = 0g.

Theorem 4 [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF], [START_REF] Zhao | On su¢ciency of the Kuhn-Tucker conditions in nondi¤erentiable programming[END_REF] (Generalized Karush-Kuhn-Tucker necessary optimality conditions). Let x 2 D be an optimal solution in problem (P) and some suitable constraint quali…cation be satis…ed at x. Then, there exist 2 R m , 2 R s such that

0 2 @f (x) + m X i=1 i @g i (x) + s X j=1 j @h j (x), (1)
i g i (x) = 0, i 2 J; (2)
i 2 R + ; i 2 J: (3)
We will assume that a suitable constraint quali…cation is satis…ed at any optimal point in problem (P).

De…nition 5 The point x 2 D is said to be Karush-Kuhn-Tucker point (a KKT point, for short) if there exist the Lagrange multipliers 2 R m , 2 R s such that the conditions (1)-( 3) are satis…ed at x.

The Exactness of the Exact l 1 Penalty Function Method

The most popular nondi¤erentiable exact penalty function is the absolute value penalty function also called the exact l 1 penalty function. Its de…nition, for the considered optimization problem (P), is the following

minimize P (x; c) = f (x) + c 2 4 X i2I g + i (x) + X j2J jh j (x)j 3 5 ; (P(c)) (4)
where, for a given constraint g i (x) 0, the function g + i is de…ned by

g + i (x) = 0 if g i (x) 0; g i (x) if g i (x) > 0: (5)
The unconstrained optimization problem de…ned above, we call the penalized optimization problem with the absolute value penalty function.

It is known (see, for example, [START_REF] Bazaraa | Nonlinear programming: theory and algorithms[END_REF]) that under suitable convexity assumptions and a constraint quali…cation, there exists a …nite value c that will recover an optimal solution in the constrained optimization problem (P) via the minimization of the exact penalty function being the objective function in the exact penalized optimization problem (P(c)). Now, we generalize this result by weakening the convexity assumption imposed on the functions constituting the considered nonsmooth optimization problem (P).

Theorem 6 Let x 2 D be a Karush-Kuhn-Tucker point in the constrained optimization problem (P), at which the Generalized Karush-Kuhn-Tucker conditions (1)-( 3) are satis…ed with the Lagrange multipliers 2 R m and 2 R s . Let J + (x) = j 2 J : j > 0 and J (x) = j 2 J : j < 0 . Furthermore, assume that the functions f , g i , i 2 I, h j , j 2 J + (x), are locally Lipschitz invex at x on X with respect to the same function and the functions h j , j 2 J (x), are locally Lipschitz incave at x on X with respect to the same function . If c is assumed to be su¢ciently large (it is su¢cient to set c max i , i 2 I, j , j 2 J , where i , i = 1; :::; m; j , j = 1; :::; s, are the Lagrange multipliers associated with the constraints g i and h j , respectively), then x is also a minimizer of its penalized optimization problem (P(c)) with the absolute value penalty function.

Proof. By assumption, x is a Karush-Kuhn-Tucker point in the constrained optimization problem (P), at which the Generalized Karush-Kuhn-Tucker conditions (1)-( 3) are satis…ed with the Lagrange multipliers 2 R m and 2 R s . Since c max i , i 2 I, j , j 2 J , then, by de…nition of the objective function in the penalized optimization problem (P(c)), it follows that

P (x; c) = f (x)+c m X i=1 g + i (x)+c s X j=1 jh j (x)j f (x)+ m X i=1 i g + i (x)+ s X j=1 j h j (x) : (6
) Thus, [START_REF] Craven | Invex functions and constrained local minima[END_REF] gives

f (x) + m X i=1 i g + i (x) + s X j=1 j h j (x) f (x) + m X i=1 i g i (x) + s X j=1 j h j (x). (7)
By assumption, the inequality constraints g i , i 2 I, and the equality constraints h j , j 2 J + (x), are locally Lipschitz invex at x on X the equality constraints h j , j 2 J (x), are locally Lipschitz incave at x on X. Hence, by the Generalized Karush-Kuhn-Tucker conditions (2) and (3) together with the feasibility of x in problem (P), it follows that the inequality

f (x)+ m X i=1 i g i (x)+ s X j=1 j h j (x) f (x)+ m X i=1 i T i (x; x)+ s X j=1 j T j (x; x) (8)
holds for any i 2 @g i (x), i = 1; :::; m, and for any j 2 @h j (x), j = 1; :::; s. Then, using the Generalized Karush-Kuhn-Tucker condition (2) , we get

f (x) + P m i=1 i h g i (x) + T i (x; x) i + P s j=1 j h j (x) + T j (x; x) = f (x) + P m i=1 i T i (x; x) + P s j=1 j T j (x; x) : (9)
Thus, by the Generalized Karush-Kuhn-Tucker necessary optimality condition (1), it follows that

f (x) + m X i=1 i T i (x; x) + s X j=1 j T j (x; x) = f (x) T (x; x) , (10)
where 2 @f (x). By assumption, f is locally Lipschitz invex at x on X also with respect to the function . Using De…nition 2 together with the feasibility of x in problem (P), we get

f (x) T (x; x) f (x) = f (x) + c m X i=1 g + i (x) + c s X j=1 jh j (x)j = P (x; c) . ( 11
)
Then, by ( 6)-( 11), we conclude that the inequality P (x; c) P (x; c) holds for all x 2 X. This means that x is a minimizer of the penalized optimization problem (P(c)) with the absolute value penalty function and the proof of theorem is complete.

Corollary 7 Let x be an optimal point in the considered optimization problem (P). Furthermore, assume that all hypotheses of Theorem 6 are ful…lled. Then x is also a minimizer in the penalized optimization problem (P(c)) with the absolute value penalty function.

Theorem 8 Let the point x be a minimizer of the penalized optimization problem (P(c)) with the absolute value penalty function . Furthermore, assume that the functions f , g i , i 2 I, h j , j 2 J + (e x), are locally Lipschitz invex at e x on X with respect to the same function , and the functions h j , j 2 J (e x), are locally Lipschitz incave at e

x on X with respect to the same function , where e

x is any Karush-Kuhn-Tucker point in problem (P), at which the Karush-Kuhn-Tucker necessary optimality conditions ( 1)-( 3) are satis…ed with the Lagrange multipliers e 2 R m and e 2 R s . If the set of all feasible solutions in the constrained optimization problem (P) is compact and the penalty parameter c is su¢ciently large (it is su¢cient if c satis…es the following condition c > max n e i , i 2 I, e j , j 2 J o ), then x is also optimal in problem (P).

Proof. We assume that x is a minimizer in the penalized optimization problem (P(c)) with the absolute value penalty function. Then, by the de…nition of the penalized optimization problem (P(c)) and ( 5), the following inequalities

f (x)+c P m i=1 g + i (x) + P s j=1 jh j (x)j f (x)+c P m i=1 g + i (x) +
P s j=1 jh j (x)j f (x) hold for all x 2 X. Thus, for all x 2 D, the following inequality

f (x) f (x) (12)
holds. The inequality above means that values of the function f are bounded below on the set D of all feasible solutions in the constrained optimization problem (P). Since f is a continuous function bounded below on the compact set D, therefore, by Weierstrass' theorem, f admits its minimum e x on D. Now, we prove that x is also optimal in the considered optimization problem (P). First, we show that x is feasible in problem (P). By means of contradiction, suppose that x is not feasible in problem (P). As we have established above, the given constrained optimization problem (P) has an optimal solution e

x. Since a constraint quali…cation is satis…ed at e

x, then there exist the Lagrange multipliers e 2 R m and e 2 R s that the Generalized Karush-Kuhn-Tucker necessary optimality conditions (1)-( 3) are satis…ed at e

x. By assumption, the functions f , g i , i 2 I, h j , j 2 J + (e x), are invex at e x on X with respect to the same function and the functions h j , j 2 J (e x), are incave at e

x on X with respect to the same function . Therefore, by De…nition 2, respectively, it follows that the inequalities f (x) f (e x) T (x; e x) ;

g i (x) g i (e x) T i (x; e x) ; i 2 I, (13)
h j (x) h j (e x) T j (x; e x) ; j 2 J + (e x), (14)
h j (x) h j (e x) T j (x; e x) ; j 2 J (e x) (15)
hold for each 2 @f (e x), i 2 @g i (e x), i = 1; :::; m, and j 2 @h j (e x), j = 1; :::; s. Multiplying ( 14), ( 15) and ( 16) by the associated Lagrange multiplier and then adding both sides of the obtained inequalities and both sides of ( 13), we get

f (x) f (e x) + P m i=1 e i g i (x) P m i=1
e i g i (e x) + P s j=1 e j h j (x)

P s j=1 e j h j (e x) h T + P m i=1 e i T i + P s j=1 e j T j i (x; e x) :
Using ( 5) with the Karush-Kuhn-Tucker necessary optimality conditions (1), ( 2) and the feasibility of e x in problem (P), we get

f (x) + m X i=1 e i g + i (x) + s X j=1 e j jh j (x)j f (e x) . (17)
By assumption, the penalty parameter c is su¢ciently large (it is su¢cient that c > max n e i , i 2 I, e j , j 2 J o . Since x is assumed to be not feasible in the given optimization problem (P), therefore, at least one of g + i (x) and jh j (x)j must be nonzero. Therefore, (17) yields

f (x) + c 2 4 m X i=1 g + i (x) + s X j=1 jh j (x)j 3 5 > f (e x) . (18)
Then, by e x 2 D and (2), we get

f (x) + c 2 4 m X i=1 g + i (x) + s X j=1 jh j (x)j 3 5 > f (e x) + c 2 4 m X i=1 g + i (e x) + s X j=1 jh j (e x)j 3 5 .
Then, by the de…nition of the exact l 1 penalty function (see ( 4)), it follows that the following inequality P (x; c) > P (e x; c) holds, which is a contradiction to the assumption that x is a minimizer in the penalized optimization problem (P(c)) with the absolute value penalty function. Thus, we have proved that x is feasible in the given constrained optimization problem (P). Hence, the optimality of x in problem (P) follows directly from [START_REF] Kim | On duality theorems for nonsmooth Lipschitz optimization problems[END_REF].

Corollary 9 Let the hypotheses of Corollary 7 and Theorem 8 are ful…lled.

Then, the set of optimal solutions in the considered extremum problem (P) and the set of minimizers in its associated exact penalized optimization problem (P(c)) with the absolute value penalty function coincide.

Example 10 Consider the following nonsmooth optimization problem

f (x) = arctan (jxj) ! min g(x) = 1 2 e jxj x 1 0: (P1)
Note that D = fx 2 R : x 0g and x = 0 is an optimal solution in the considered nonsmooth optimization problem (P1). Since we use the exact l 1 penalty method for solving problem (P1), then we construct the following unconstrained optimization problem P (x; c) = arctan (jxj) + c max 0; 1 2 e jxj x 1 ! min : (P1(c))

Note that x = 0 is feasible in problem (P1) and the Generalized Karush-Kuhn-Tucker necessary optimality conditions (1)-( 3) are at x with the Lagrange multiplier satisfying the following condition: 0 2 @f (x) + @g (x), where @f (x) = [ 1; 1] and @g (x) = [ 1; 0]. Further, it can be established by De…nition 2 that the objective function f and the constraint function g are locally Lipschitz invex at x on R with respect to the same function de…ned by (x; x) = 1 2 (arctan (jxj) arctan (jxj)). Then, by Theorems 6 and 8, it follows that, for any penalty parameter c satisfying c > , there is the equivalence between the sets of optimal solutions in optimization problems (P1) and (P1(c)). Further, note that not all functions involved in problem (P1) are di¤erentiable and convex. Therefore, in order to show that the point x = 0, being optimal in (P1), is also a minimizer in the unconstrained optimization problem (P1(c)), we can not use the conditions for convex smooth optimization problem (see, for instance, Theorem 9.3.1 [START_REF] Bazaraa | Nonlinear programming: theory and algorithms[END_REF]).

Example 11 Consider the following nonsmooth constrained optimization problem

min f (x) = 8 > > < > > : x + 4 if x < 4; 1 2 x + 10 if 4 x < 0; 5x + 10 if 0 x < 2; x 2 if x 2; g(x) = x 1 4 0, (P2)
in which not all functions are locally Lipschitz invex. Note that D = x 2 R : x 1 4 and x = 4 is an optimal solution in the considered optimization problem (P2). Since 0 2 @f (0) = 5; 1 2 , then e x = 0 is a stationary point of f . It is not di¢cult to show that e

x is not a global minimizer of f . Then the objective function f is not locally Lipschitz invex on R with respect to any function de…ned by : R R ! R (see, for example, [START_REF] Reiland | Nonsmooth invexity[END_REF]). However, we use the exact l 1 penalty method to solve the considered optimization problem (P2). Therefore, we construct the following unconstrained optimization problem

P (x; c) = f (x) + c max 0; x 1 4 ! min (P2(c))
Note that x = 4, being an optimal solution in problem (P2), is not a global minimizer in the associated penalized optimization problem (P2(c)) for all values of the penalty parameter c satisfying the condition c > = 0, (where is the Lagrange multiplier associated with the inequality constraint g satisfying the Karush-Kuhn-Tucker necessary optimality conditions (1)-( 3)). However, for every penalty parameter c 2 (0; 32 7 ), the point b x = 2 is a global minimizer in the above penalized optimization problem (P2(c)). Therefore, there is no the equivalence between the sets of optimal solutions in problems (P2) and (P2(c)) for any penalty parameter c satisfying the condition c > . This follows from the fact that not all functions constituting the considered optimization problem (P2) are locally Lipschitz invex on R.

Remark 12 Peressini et al. [START_REF] Peressini | The mathematics of nonlinear programming[END_REF] considered di¤erentiable convex optimization problems and solved them by using the exact l 1 penalty function method. Under assumption that the objective function in the constrained optimization problem is coercive (see De…nition 3), they proved that, for su¢ciently large values of the penalty parameter c, the constrained optimal solution in (P) is also a minimizer in its associated penalized optimization problem (P(c)) with the exact l 1 penalty function. But the …nite value of the penalty parameter c, above which this result holds, was not given in [START_REF] Peressini | The mathematics of nonlinear programming[END_REF]. Note that the objective function in the optimization problem (P2) considered in Example 11 is coercive. However, for not all values of the penalty parameter c satisfying the condition c > , an optimal solution in the considered optimization problem (P2) yields a minimizer in its associated penalized optimization problem (P(c)) with the exact l 1 penalty function. But the result proved in the paper shows that, under invexity assumptions imposed on the functions constituting the constrained nonsmooth optimization problem (P), for every value of the penalty parameter c satisfying the condition c > max i , i 2 I, j , j 2 J , the sets of optimal solutions in problems (P) and (P(c)) coincide. Hence, this example shows that in the case when the objective function is coercive but not invex the result established in the paper might not be true for such optimization problems.


Marek D Lugosz
email: mdlugosz@agh.edu.pl

The Minimum Energy Building Temperature Control

Keywords: building temperature model, control, optimization

One of the most important factors of users comfort inside building is air temperature. From the other side one of the most biggest position in home budget is price for heat. Mutually exclusive indices are the cause that the control of temperature task using the smallest amount of energy as it is possible is very difficult. In this paper is presented simple model of temperature changes inside building base on lumped capacity method. Using this method finally obtains mathematical model of temperature changes which model is equivalent in structure to electrical RCnetwork. The model is composed of linear differential equations. Based on this mathematical model the simple algorithm controlling of temperature inside room is proposed. In this article are also included numerical simulations of the proposed solutions.

Introduction

Problem of optimal use of the heat energy for heating residential building is still a current problem. Some of the main reasons for this are: still rising energy price (electricity, gas, coal), still rising power consumption by household or environmental pollution. The goal of this article is to present a control system which stabilises the temperature inside the building with the using minimum amount of energy.

The air temperature inside the building T i depends on many factors. Some of them like: solar radiation, wind, heating system, light, people, air ventilation are showed on figure 1. Some of these factors are unpredicted like: people inside, light, air ventilation. Some of them are periodical and can be measured or predicted, for example: solar radiation, temperature outside, wind direction and force. The physical phenomena of thermal conductivity are also very complex and described by partial differential equations which depend on time and spatial variables. For those reasons one and general thermal model of the building does not exist. On the other hand for searching optimal controls the mathematical model of the Fig. 1. The temperature inside the building depends on many factors such as: heating system, external air temperature, wind, solar heat, casual heat gains, structure of the building.

system is necessary. In the literature can be found three main methods to obtain and identify an approximated thermal model of the building.

The first method: the impulse response factor method [START_REF] Clarke | Energy simulation in building design[END_REF][START_REF] Underwood | Modelling methods for energy in buildings[END_REF] is based on the response of the model if the excitation is a unit impulse. Making some additional assumptions and using the properties of Laplace transform the response of the wall to this excitation function can be expressed as time series.

The second method is the finite difference method. This is a numerical method for solving partial differential equation of the heat conduction [START_REF] Clarke | Energy simulation in building design[END_REF][START_REF] Underwood | Modelling methods for energy in buildings[END_REF]. The finite difference method is based on approximation derivates by algebraic equation. The building wall is divided into a finite number of layers and temperature for each layers is computed using set of the algebraic equations.

The third method: the lumped parameter method (or other name the lumped capacitance method) base on assumptions that transfer of the heat flux between two spaces which are divided by partition (wall) can be modeled by the equivalent electrical RC circuit [START_REF] Gouda | Low-order model for the simulation of a building and its heating system[END_REF][START_REF] Gouda | Building thermal model reduction using nonlinear constrained optimization[END_REF][START_REF] Underwood | Modelling methods for energy in buildings[END_REF]. The parameters of the electrical RC circuit like resistances are interpreted as thermal resistances, capacities are interpreted as heat capacities of the modeled elements. The physical properties of the construction elements of the building are represented by resistors and capacitors. The lumped parameter method describes changes of the air temperatures or the temperatures of the construction elements in one point. Finally, the mathematical model which is obtained by using the lumped parameter method has the form of linear differential equations. This model can be easily solved by analytical or numerical methods. In this paper, the lumped parameter method was chosen for modeling changes of the indoor air temperature of building.

The plan of the article is as follow, the first section contains short description of the lumped parameter method (LPM), next section contains description of the LQR controller. The last section presents some of experimental results. At the end of the article are contained conclusions and plans for the future works.

Thermal modeling methodologies

The most suitable form of the mathematical model of dynamic system for searching of optimal control solutions, is form of the linear differential equations.

ẋ(t) = Ax(t) + Bu(t) + Zz(t) (1)
y(t) = Cx(t) (2)
where A n×n -state space matrix, B n×r -control matrix, Z n×k -noise matrix,

C m×n -output matrix, x(t) ∈ X = R n -state space vector, z(t) ∈ Z = R k - noise vector, u(t) ∈ U = R r -control vector, y(t) ∈ Y = R m -output vector.
In all simulations the noise matrix Z was assumed to zero.

LPM

The assumptions of the lumped parameters method is that the temperature of the solid is spatially uniform at any instant during the transport of the heat process [START_REF] Incropera | Fundamentals of Heat and Mass Transfer[END_REF]. The result of this assumption is that the heat flow between two spaces which are separated by partition can be replaced by an equivalent RC electrical circuit [START_REF] Gouda | Low-order model for the simulation of a building and its heating system[END_REF][START_REF] Gouda | Building thermal model reduction using nonlinear constrained optimization[END_REF][START_REF] Underwood | Modelling methods for energy in buildings[END_REF]. The lumped parameter method describes changes of the temperature in one point so it is only an approximation of the real temperature. These simplifications allow us use the linear differential equations instead of more complicated partial differential equations. The lumped parameters method can be used for materials for which the conductivity in the middle is larger than the conductivity on the material surface [START_REF] Crabb | A simplified thermal response model[END_REF][START_REF] Incropera | Fundamentals of Heat and Mass Transfer[END_REF].

As was said, the heat flow between two spaces which are separated by the partition can be replaced by the electrical circuit and figure 2 shows this. The meaning parameters are: the R out and R int thermal resistances of area outer and inner, C total thermal capacity of the partition. The equation of the heat conduction based on the first-order model is:

C total dT dt = (T o -T ) R out + (T i -T ) R in + q (3)
where q represents the other heat sources, T is uniform material temperature, T o is outer air temperature and T i is inner air temperature. The wall on figure 2 consists one of the uniform material but in the real wall may be build more than one of layers of the uniform materials. In this case, we can extend the model by adding the next equations for each uniform layer [START_REF] Gouda | Building thermal model reduction using nonlinear constrained optimization[END_REF][START_REF] Underwood | Modelling methods for energy in buildings[END_REF]. Also, if is needed the more accurate mathematical model we may add the next equations to the model [START_REF] Deng | Building thermal model reduction via aggregation of states[END_REF][START_REF] Gouda | Building thermal model reduction using nonlinear constrained optimization[END_REF]. All these operations finally increase total number of the equations and order of the model.

Ti

The simple thermal model of building

The state equations from ( 4) to (9) describe the thermal behaviour of more complex space [START_REF] Gouda | Low-order model for the simulation of a building and its heating system[END_REF]. The space inside which temperature is modeled contains two external walls , two partitions, floor and ceiling. In this case changes of the temperature of the indoor air depend on much more factors.

C 1 Ṫ1 = U 1 (T i -T 1 ) + U 2 (T 0 -T 1 ) (4)
C 2 Ṫ2 = U 3 (T i -T 2 ) + U 4 (T 0 -T 2 ) (5)
C 3 Ṫ3 = U 5 (T i -T 3 ) + U 6 (T z1 -T 3 ) + Q s (6) C 4 Ṫ4 = U 6 (T i -T 4 ) + U 5 (T z2 -T 4 ) (7)
C 5 Ṫ5 = U 7 (T i -T 5 ) + U 7 (T z3 -T 5 ) (8)
C 6 Ṫi = U 1 (T 1 -T i ) + U 3 (T 2 -T i ) + U 5 (T 3 -T i ) (9) +U 6 (T 4 -T i ) + U 7 (T 5 -T i ) + U 8 (T o -T i ) +Q p + Q g
The parameters are: T 1 and T 2 -temperature of the building structure, T 3 and T 4 -temperature of the floor and ceiling, T 5 -temperature of the partitions, T o -outdoor air temperature and T i -indoor air temperature. The figure 3 shows electrical circuit RC which is equivalent with the building thermal model. The electrical parameters of this circuit correspond with physical parameters of the building. The resistances are equivalent to overall thermal transmittance and capacities are equivalent to thermal capacity. As is shown on the figure 3 electrical circuit has the form of the RC ladder network. The analytical solutions of the Fig. 3. The equivalent electrical circuit RC for the building thermal model given by equations number from (4) to (9) model's equations can be easily found and analysed and this is big advantage of this type of the model.

In next simulations values of the parameters of the building like C i and R i are the same as those adopted in the article [START_REF] Gouda | Low-order model for the simulation of a building and its heating system[END_REF].

Temperature control in the building

The typical idealised behaviour of indoor air temperature is shown on figure 4. As we can see the three different phases can be highlighted. The first phase, when indoor air temperature should reach the reference value in given time. The second phase when the air temperature should be stabilised on the specified level. The third phase when the air temperature do not need to be stabilised or controlled. This is idealised behaviour of indoor air temperature but generally all more complicated schemas of the temperature changes can be described by using those three phases.

Only in the first and second phases is required active control of indoor air temperature. The main goal of the control system is to control indoor temperature in the first and second phases but in the first phase, time of control is also limited. In this paper is proposed to use two different controllers. The finitehorizon LQR controller which works in the first phase and the infinite-horizon LQR controller which works in the second phase. The next two subsections are describe shortly those controllers and are present their advantages and disadvantages.

The finite-horizon LQR controller

The finite-horizon LQR controller minimizes the cost function (10) [START_REF] Mitkowski | Stabilizacja Systemów Dynamicznych[END_REF]:

J(x, u) = 1 2 t k 0 x(t) T W x(t) + u(t) T Ru(t) dt + 1 2 x(t k ) T F x(t k ). ( 10
)
The matrices W , F , R are weight matrices and that matrices must be nonnegative and symmetric and W = W T 0, F = F T ≥ 0 R = R T > 0, the pair of matrices (A, B) is stabilisabe , the pair of matrices (W , A) is detectable, t kis the control time. The control law is given by equation (11) [START_REF] Mitkowski | Stabilizacja Systemów Dynamicznych[END_REF]:

u(t) = -R -1 B T K(t)x(t) (11)
where matrix K is unique, symmetric and nonnegative solution of Riccati differential equation (12) [START_REF] Mitkowski | Stabilizacja Systemów Dynamicznych[END_REF]:

K(t) = K(t)BR -1 B T K(t) -A T K(t) -K(t)A -W . ( 12
)
The controller (11) is nonstationary because the values of matrix K depend on time.

Simulation The fig. 5 shows the result of simulation of the control system with the finite-horizon LQR controller. The control task was to increase the temperature value from 15 to 20 degrees in a finite time (75 minutes). The first plot shows the change of the indoor air temperature T i (9), the second graph shows the change of the control signal u(t) (11). The finite-horizon LQR controller is complicated in practical applications. First of all this is the nonstationary controller because the gain matrix K(t) is depends on time. In order to compute matrix K(t) the Riccati differential equation ( 12) must be solved. The same result, raise the value of indoor temperature from one level to other level in finite time, can be obtained using the infinite-horizon LQR controller with appropriate chose of the weight matrices, W and K.

The infinite-horizon LQR controller

The infinite-horizon LQR controller minimizes a cost function (13) [START_REF] Mitkowski | Stabilizacja Systemów Dynamicznych[END_REF]:

J(x, u) = ∞ 0 x(t) T W x(t) + u(t) T Ru(t) dt (13)
where W = W T ≥ 0, R = R T > 0, the pair of matrices (A, B) is stabilisable, the pair of matrices (W , A) is detectable. The cost function (13) contains two parts: a part which is joined with state space vector x(t) and part which is joined with the control vector u(t). The matrices W and R are called the weight matrices and appropriate selection of their values determines which part of cost function is better stabilised. The infinite-horizon LQR controller is proportional controller and the control law is given by equation (14) [START_REF] Mitkowski | Stabilizacja Systemów Dynamicznych[END_REF]:

u(t) = -R -1 B T Kx(t) (14)
where K is unique, symetric, nonnegative solution of algebraic Riccati equation:

KBR -1 B T K -A T K -KA -W = 0. (15)
Simulation The fig. 6 shows the result of simulation of the control system with the infinite-horizon LQR controller. The control task was to increase the temperature value from 15 to 20 degrees in a finite time (1 hour). The first plot of figure 6 shows the change of the indoor temperature T i (9), the second graph shows the change of the feedback control signal u(t) (11). The infinite-horizon LQR controller is easier in practical applications because the values of matrix K are constant and it is a stationary proportional controller. By changing values of the matrices W and R is possible to modify in wide range of how the controller works e.g. approximate time after which the desired value of controlled variable will be achieved. 4 The complex control system

The goal is to build a control system which works properly in phase I and phase II, see fig. 4. The control system will be implemented in a computer so this gives ability to build more complex control system. The infinite-horizon LQR controller will be used because, as was said earlier, this kind of controller is easy to use in practical applications. The figure 7 shows the block diagram of proposed control system. The LQR controller number 1 works during phase I The fig. 9 shows a comparison of the controls signals for the single LQR controller (dotted line) and the control system which contains the two LQR controllers (solid line). The control system which includes two LQR controllers uses less energy to stabilise the temperature than system with one LQR controller and the area between two curves corresponds to the amount of saved energy.

Conclusions

This paper is presented the control system which contains two infinite-horizon LQR controllers. The infinite-horizon LQR controllers were chosen because: the control system is closed-loop system with negative feedback, simple structure of the infinite-horizon LQR controller (proportional controller) and modifying the weight matrices W and R can change the nature of the work control system. Also, some disadvantages of the infinite-horizon LQR controller are existing like: the LQR controller is a proportional controller so always is a deviation between desired value and real value of controlled signal, the LQR controller for proper work needs to know values of all coordinates of the state vector x(t) and in some cases reconstruction of the non-measurable coordinates of the state variables is

Fig. 2 .

 2 Fig. 2. Representation of the lumped parameter construction element

Fig. 4 .

 4 Fig. 4. Typical idealised changes of indoor air temperature.

Fig. 5 .

 5 Fig. 5. (a): Changes of the indoor air temperature Ti (9) and (b): the control signal u(t) = Qp(t) (11).

Fig. 6 .

 6 Fig. 6. (a): Changes of the inside temperature Ti (9) and (b): the control signal u(t) = Qp(t) (11).

Fig. 7 .

 7 Fig. 7. Scheme of the proposed complex control system

Fig. 8 .

 8 Fig. 8. Changes of the inside temperature (left plot) and the control signal u(t) = Qp(t) (right plot).

Fig. 9 .

 9 Fig. 9. Comparison of the controls signals u(t) = Qp(t) for the single LQR controller (dotted line) and the control system which contains the two LQR controllers (solid line).

  needed. As shown by results of the simulations there is possibility to control indoor air temperature efficiently and using less energy. The LQR controller minimises the cost function which also takes into account the energy consumption of the control signal. Recent times, can be observed the growing popularity of wireless home automation devices. The future work will be concentrated on practical implementation of the proposed solutions in devices which work in ZWave and ZigBee standard.Work is financed by NCN-National Science Centre founds for 2011-2013 as a research project. Contract number N N514 644440.


F.-Javier Heredia
email: f.javier.heredia@upc.edu

Cristina Corchero
email: cristina.corchero@upc.edu

Eugenio Mijangos
email: eugenio.mijangos@ehu.es

Solving electric market quadratic problems by Branch and Fix Coordination methods ⋆

Keywords: Liberalized Electricity Market, Optimal Bid, Stochastic Programming, Quadratic Branch-and-Fix Coordination

The electric market regulation in Spain (MIBEL) establishes the rules for bilateral and futures contracts in the day-ahead optimal bid problem. Our model allows a price-taker generation company to decide the unit commitment of the thermal units, the economic dispatch of the bilateral and futures contracts between the thermal units and the optimal sale bids for the thermal units observing the MIBEL regulation. The uncertainty of the spot prices is represented through scenario sets. We solve this model on the framework of the Branch and Fix Coordination metodology as a quadratic two-stage stochastic problem. In order to gain computational efficiency, we use scenario clusters and propose to use perspective cuts. Numerical results are reported.

Introduction

This work is applied to the Iberian Electricity Market (MIBEL) comprising the Spanish and Portuguese electricity systems. The MIBEL market includes in the short-term: the day-ahead market (DAM) and a set of balancing, reserve and adjustment markets (intraday markets); these markets are complemented with the medium-and long-term mechanisms: a derivatives market and different kinds of bilateral contracts. This structure is similar to other European electricity markets and explains why generation companies can no longer optimize their short-term generation planning decisions, i.e. their bidding strategies, without considering the relationship between the short-term bid and the medium-term physical products. The MIBEL's directives dictate specific rules describing how these medium-term mechanisms should be included into the DAM bid. This work deals with the most relevant medium-term mechanisms in the MIBEL: the national bilateral contracts (BC) and the future physical contracts (FC). Stochastic programming techniques are applied to maximize the expected value of the utility's profit coming from the day-ahead market, where the significative random variable is the auction clearing price of the day-ahead electricity market. This random variable is modeled through a set of scenarios of the forecasted prices. The set of scenarios is used to feed a two-stage stochastic optimization model that finds the optimal day-ahead bid of a price-taker GenCo (an electrical utility without influence over the market prices) operating in the MIBEL and holding bilateral and physical futures contracts.

The extensive form of the deterministic equivalent of this stochastic programming problem will be a mixed integer quadratic programming problem (MIQP), which is difficult to solve efficiently, particularly for large-scale instances. Several algorithmic approaches can be adopted to overcome this difficulty. In [START_REF] Corchero | A new optimal electricity market bid model solved through perspective cuts[END_REF] the quadratic objective function of this problem is approximated by a polyhedral outer approximation by means of perspective cuts so that we can exploit the efficiency of general-purpose solvers for mixed integer linear problems (MILP). An alternative to the perspective cuts methodology is the Second-Order Cone Program reformulation (SOCP, [START_REF] Tawarmalani | Semidefinite relaxations of fractional programs via novel convexification techniques[END_REF]), but for quadratic problems the perspective cuts reformulation was reported to be more efficient [START_REF] Frangioni | A computational comparison of reformulations of the perspective relaxation: SOCP vs. cutting planes[END_REF]. Finally, the Branch-and-Fix Coordination (BFC) method has been used successfully to solve two-stage stochastic mixed integer linear problems [START_REF] Escudero | A general algorithm for solving two-stage stochastic mixed 0-1 first-stage problems[END_REF] to solve the day-ahead optimal bid problem. In this work we propose an combination between BFC and PC to efficiently solve the optimal day-ahead bid problem.

2 Day-ahead electricity market bid with futures and bilateral contracts model (DAMB-FBC)

In this section the model (DAMB-FBC) is formulated as a two-stage stochastic programming problem that allows a price-taker generation company to optimally decide the unit commitment of its thermal units, the economic dispatch of the bilateral and futures contracts between the thermal units, and the optimal generation bid of the committed units to the MIBEL's day-ahead market. The objective function of the model represents the expected profits of the GenCo's participation in the day-Ahead market. The constraints assure that the MIBEL's rules and the operational restrictions of the units are respected. The main decision variables are the ones that model the start-up and shut-down of the units, the quantity that will be bid at instrumental price and the scheduled energy committed to the bilateral and the futures contracts settlement.

Parameters

The (DAMB-FBC) model considers a price-taker GenCo owning a set of thermal generation units I that bid to the t ∈ T = {1, 2, .., 24} hourly auctions of the DAM. The parameters for the i th thermal unit are:

c b i , c l i and c q i , generation costs with constant, linear and quadratic coefficients (e, e/MWh and e/MWh 2 respectively).

-P i and P i , upper and lower bounds on the hourly energy generation (MWh).

c on i and c of f i , start-up and shut-down costs (e). -t on i and t of f i , minimum operation and minimum idle time (h).

A base load physical futures contract j ∈ F is defined by:

-U j , the set of generation units allowed to cover the FC j.

-L F j , the amount of energy (MWh) to be procured each interval of the delivery period by the set U j of generation units to cover contract j.

λ F j , the price of contract j (e/MWh).

A base load bilateral contract k ∈ B is defined by:

-L B k , the amount of energy (MWh) to be procured at each interval of the delivery period by the set of available generation units to cover the BCs.

λ B k , the price of the contract k (e/MWh).

The random variable λ D t , the clearing price of the t th hourly auction of the DAM, is represented in the two-stage stochastic model by a set of scenarios s ∈ S, each one with its associated clearing price for each DAM auction t ∈ T :

λ D,s t clearing price for auction t at scenario s (e/MWh).

-P s probability of scenario s.

Variables

Those decision variables that doesn't depend on the scenarios are called first stage (or here-and-now) variables and in our formulation are, for each t ∈ T and i ∈ I:

u ti , the unit commitment (binary) -c u ti , c d ti , the start-up/shut-down costs variables. -q ti , the instrumental price offer bid.

f tij , the scheduled energy for FC j ∈ F.

b ti , the scheduled energy for the pool of BCs .

Decision variables that can adopt different values depending on the scenario are called second stage variables and in our formulation are, for each t ∈ T , i ∈ I and scenario s ∈ S:

g s ti , the total generation. -p s ti , the matched energy in the day-ahead market.

Contraints

Bilateral and futures contracts constraints The coverage of both the physical futures contracts and the bilateral contracts must be guaranteed. The constraints for each futures contract are:

i∈Uj f tij = L F j t ∈ T , j ∈ F (1)
f tij ≥ 0 t ∈ T , j ∈ F, i ∈ I (2)
and the bilateral contract constraints are:

i∈I b ti = k∈B L B k t ∈ T (3) 0 ≤ b ti ≤ P i u ti i ∈ I, t ∈ T (4)
where L B k is the energy that has to be settled for contract k ∈ B Day-ahead market and total generation constraints As we have introduced, we will use the value of the matched energy in our formulation. The matched energy is the accepted energy in the clearing process, that is, the energy generated that will be rewarded at the clearing price. This matched energy is uniquely determined by the sale bid and the clearing price and it will play a central role in the presented model [START_REF] Corchero | A new optimal electricity market bid model solved through perspective cuts[END_REF]. The MIBEL's rules affecting the day-ahead market establishes the relation between the variables representing the matched energy p s ti , the energy of the bilateral contracts b ti , the energy of the futures contracts f tij , the instrumental price offer bid q ti , and the commitment binary variables u ti . The energies L F j and L B k must be integrated in the MIBEL's DAM bid observing the two following rules:

1. If generator i contributes with f tij MWh at period t to the coverage of the FC j, then the energy f tij must be offered to the pool for free (instrumental price bid). 2. If generator i contributes with b ti MWh at period t to the coverage of any of the BCs, then the remaining production capacity P ib ti must be bid to the DAM.

These rules can be included in the model by means of the following set of constraints:

p s ti ≥ q ti i ∈ I, t ∈ T , s ∈ S (5)
p s ti ≤ P i u ti -b ti i ∈ I, t ∈ T , s ∈ S (6)
q ti ≥ P i u ti -b ti i ∈ I, t ∈ T , s ∈ S (7)
q ti ≥ j | i∈Uj f tij i ∈ I, t ∈ T , s ∈ S (8)
where:

(5) and ( 6) ensure respectively that the matched energy p s ti will be greater than the instrumental price bid q ti and less than the total available energy not allocated to BC. ( 7) and ( 8) guarantee respectively that the instrumental price bid will be greater than the minimum generation output of the unit and greater than the contribution of the unit to the FC coverage.

Please note that (2) together with [START_REF] Mijangos | An algorithm for two-stage stochastic quadratic problems[END_REF] assures that q ti will be always nonnegative. The total generation level of a given unit i, g s ti , is defined as the addition of the allocated energy to the BC plus the matched energy of the DAM:

g s ti = b ti + p s ti i ∈ I, t ∈ T , s ∈ S (9)
Contraints ( 1)-( 9) assure that g s ti will be always either zero or g s ti ∈ [P i , P i ], that is:

P i u ti ≤ g s ti ≤ P i u ti , i ∈ I, t ∈ T , s ∈ S (10)
Unit commitment constraints This section includes the formulation for the unit commitment of the thermal units [START_REF] Corchero | A new optimal electricity market bid model solved through perspective cuts[END_REF]. The first two sets of constraints model the start-up and shut-down costs and the next ones control minimum operation and idle time for each unit. First, the start-up and shut-down costs are modeled:

c u ti ≥ c on i [u ti -u (t-1)i ] i ∈ I, t ∈ T \ {1} (11)
c d ti ≥ c of f i [u (t-1)i -u ti ] i ∈ I, t ∈ T \ {1} (12)
c u ti , c d ti ≥ 0 i ∈ I, t ∈ T (13) u ti ∈ {0, 1} i ∈ I, t ∈ T (14)
The initial state of each thermal unit i can be taken into account through the parameters G i and H i that represent, respectively, the number of the initial time periods along which the thermal unit must remain on (G i ) or off (H i ). These conditions are imposed by the following set of constraints:

t+t on i -1 n=t u in ≥ t on i [u ti -u (t-1)i ] t = G i + 1, . . . , |T | -t on i + 1, i ∈ I (16) t+t of f i -1 n=t (1 -u ni ) ≥ t of f i [u (t-1)i -u ti ] t = H i + 1, . . . , |T | -t of f i + 1 i ∈ I (17)
and for the last t on i -1 and t of f i -1 time periods:

|T | n=t (u ni -[u ti -u (t-1)i ]) ≥ 0 t = |T | -t on i + 2, . . . , |T |, i ∈ I (18) |T | n=t (1 -u ni -[u (t-1)i -u ti ]) ≥ 0 t = |T | -t of f i + 2, . . . , |T |, i ∈ I (19)

Objective function

The quadratic function that gives the long-run expected profits of the GenCo after the participation in the DAM is:

min E λ D C(u, c u , c d , g, p; λ D ) = t∈T i∈I c u ti + c d ti + c b i u ti + (20) + s∈S P s (c l i g s ti + c q i (g s ti ) 2 ) -λ D,s t p s ti , (21)
where the right hand side of (20) is the on/off fixed cost of the unit commitment of the thermal units, deterministic and independent of the realization of the random variable λ D,s t and (21) represents the expected value of the benefits from the DAM. The term between parenthesis corresponds to the expression of the quadratic generation costs associated to the total generation of the unit g s ti while the last term, λ D,s t p s ti computes the incomes from the DAM due to a value p s ti of the matched energy. Please note that the constant incomes from the BC and FC, i.e. (λ F C j -λD t )L F C j , have been dropped from the objective function.

Model (DAMB-FBC)

The model defined so far can be represented as:

(DAMB-FBC)

             min E λ D C(u, c u , c d , g, p ; λ D ) s.t.
Eq. ( 1) -(4) BC and FC constraints Eq. ( 5) - [START_REF] Tawarmalani | Semidefinite relaxations of fractional programs via novel convexification techniques[END_REF] DAM and total gen. constraints Eq. ( 11)

-(19)
Unit commitment constraints Model (DAMB-FBC) is the optimization problem associated with the twostage stochastic programming problem with a set S of scenarios for the spot price λ D t , where t ∈ T . This optimization problem is a convex MIQP with a well defined global optimal solution.

QBFC method

Model (DAMB-FBC) can be rewritten as the so-called Deterministic Equivalent Model (DEM) minimize c t δ + s∈S P s q s (x, y s ) subject to :

l a ≤ A δ x ≤ u a , l s t ≤ T s   δ x y s   ≤ u s t , s ∈ S, x ≥ 0, y ≤ y s ≤ y, s ∈ S, δ ∈ {0, 1} n δ ,
where δ = u, x = (c u , c d ), y = (g, p), q(x, y) = b t x x + b t y y+ +y t Q yy y, and Q yy being a diagonal matrix.

As is showed by [START_REF] Escudero | A general algorithm for solving two-stage stochastic mixed 0-1 first-stage problems[END_REF] the compact representation (DEM) can be written as a splitting variable representation; i.e., δ and x are respectively replaced by δ s and x s , for s ∈ S. So, we have

(MIQP) minimize s∈S P s (c t δ s + q s (x s , y s )) subject to : l a ≤ A δ s x s ≤ u a , s ∈ S,
l s t ≤ T s   δ s x s y s   ≤ u s t , s ∈ S,
x s ≥ 0, y ≤ y s ≤ y, δ s ∈ {0, 1} n δ , s ∈ S,

(NAC δ ) δ s -δ s ′ = 0, ∀s, s ′ ∈ S : s = s ′ , (NAC x ) x s -x s ′ = 0, ∀s, s ′ ∈ S : s = s ′ ,
where NAC δ and NAC x are the nonanticipativity constraints.

In this method (DEM) is solved by using a Branch-and-Fix-Coordination scheme (BFC) for each scenario s ∈ S to fulfill the integrality condition (IC) on the variables δ, so that the NAC δ are also satisfied when selecting branching nodes and branching variables by the Twin-Node-Families concept (TNF), which was introduced by [START_REF] Alonso-Ayuso | BFC, a branch-and-fix coordination algorithm framework for solving some types of stochasticpure and mixed 0-1 programs[END_REF].

A similar approach to that suggested in [START_REF] Escudero | A general algorithm for solving two-stage stochastic mixed 0-1 first-stage problems[END_REF] is used in this work to coordinate the selection of the branching node and branching variable for each scenariorelated BF tree, such that the NAC δ are satisfied when fixing δ s , for all s ∈ S, either to 1 or to 0. A TNF integer set is a set of integer BF nodes (i.e. they verify IC), one per BF tree, in which the NAC δ are verified. More details about this metodology can be found in [START_REF] Mijangos | An algorithm for two-stage stochastic quadratic problems[END_REF].

When the number of scenarios is very hight, in order to gain computational efficiency we can take scenario clusters; i.e., instead a submodel for each scenario s ∈ S we can use a submodel (MIQP p ) for each scenario cluster S p ⊂ S with p = 1, . . . , p, where S p ∩ S p ′ = ∅, for all p = p ′ , and ∪ p p=1 S p = S, (MIQP p ) minimize

s∈S p P s (c t δ p + q s (x p , y s )), (22)
subject to :

l a ≤ A δ p x p ≤ u a , (23)
l s t ≤ T s   δ p x p y s   ≤ u s t , s ∈ S p , (24)
x p ≥ 0, y ≤ y s ≤ y, s ∈ S p , δ p ∈ {0, 1} n δ , (25)

These submodels are linked by the NACs δ pδ p ′ = 0 and x px p ′ = 0, for all p, p ′ ∈ {1, . . . , p} such that p = p ′ . In order to gain computational efficiency we propose to use perspective cuts (PC) [START_REF] Frangioni | Perspective cuts for a class of convex 0-1 mixed integer programs[END_REF][START_REF] Corchero | A new optimal electricity market bid model solved through perspective cuts[END_REF] to solve the quadratic subproblems in each node of the TNF. Then MIQP p becomes:

min s∈S p P s b t x x + x t Q xx x + n i=1 v s i s.t.: v s i ≥ (2q s ii y i + b s i )y s i + (c i -q s ii y 2 i )δ s i , i ∈ {1, . . . , n}, s ∈ S p , v s i ≥ (2q s ii y i + b s i )y s i + (c i -q s ii y 2 i )δ s i , i ∈ {1, .
. . , n}, s ∈ S p , Eq. ( 23) -(25).

These methods have been implemented in C++ with the help of Cplex 12.1 to solve only the quadratic subproblems. In this work two algorithmic alternatives have been considered: ⊲ QBFC: coordination of δ in the TNF of the BF trees for clusters p ∈ {1, . . . , p} without using PCs.

⊲ QBFC-PC: coordination of δ in the TNF of the BF trees for clusters p ∈ {1, . . . , p} using PCs.

For our instances the number of scenarios in each cluster is the same, |S p | = |S|/ p. Each cluster contains |S p | consecutive scenarios, starting from the first one and following in natural order.

Numerical Tests

These instances are based on the liberalized electricity market model suggested in [START_REF] Corchero | A new optimal electricity market bid model solved through perspective cuts[END_REF]. In these problems Q xx is the zero matrix, as a result, when we use perspective cuts the subproblem to solve in each node is linear. The tests have been performed on HP with Intel(R) Core(TM)2 Quad CPU Q8300 2.50GHz 4 CPU under SUSE Linux Enterprise Desktop 11 (x86 64).

In Table 1 |S| means the number of scenarios, |T | the number of periods, "# var" the number of continuous variables, "# var P CF " the number of continuous variables for the PC formulation, "# bin" the number of binary variables, and "# constr" the number of constraints for (DEM). In Table 2 below the headings QBFC are the times in CPU-seconds used for solving problems with the number of scenario cluster given below the heading p and by solving the quadratic subproblem QP p for each node using Cplex. Column QBFC-PC gives us the CPU-seconds and indicates that the quadratic subproblems QP p have been solved by using perspective cuts, which means that instead of solving a quadratic problem QP p in each node of a TNF for p ∈ {1, 2, . . . , p}, a linear problem is solved. Also, "ratio"= QBFC-PC QBFC gives us the ratio of CPU-times. Note that the running time with PC is a 30% of the running time without PC (average). The last column, "# PC", means the number of perspective cuts generated in each test.

Conclusions

We have presented an Optimal Bidding Model for a price-taker generation company operating both in the MIBEL Derivatives and Day-Ahead Electricity Market (DAMB-FBC). The model developed finds the optimal bid for the spot market, the optimal allocation of the physical futures and bilateral contracts among the thermal units and the unit commitment following in detail the MIBEL rules. The (DAMB-FBC) has been solved both with the standard BFC method and with a PC variation which reduces the running time to a 30% on the average.

( 1 -

 1 u ji ) = 0 and Hi j=1 u ji = 0, i ∈ I (15) Finally, the minimum up and down time, t on i and t of f i are imposed, up to the periods |T |-(t on i -1) and |T |-(t of f i -1), through the following set of constraints:

For

  every problem |F| = |B| = 2 and |I| = 4. If we use the PC formulation, the problem increases the number of variables in m = |T | • |I| • |S| and the number of constraints in 2 • m.

Table 1 .

 1 Test problems

	Prob. |S| |T | # var # varP CF # bin # constr
	P01 10 12 1296	1776	48	1788
	P02 20 12 2256	3216	48	3228
	P03 30 12 3216	4656	48	4668
	P04 40 12 4176	6096	48	6108
	P05 50 12 5136	7536	48	7548
	P11 10 24 2592	3552	96	3600
	P12 20 24 4512	6432	96	6480
	P13 30 24 6432	9312	96	9360
	P14 40 24 8352	12192	96	12240
	P15 50 24 10272	15072	96	15120

Table 2 .

 2 Computational results: CPU-times

	Prob. p QBFC QBFC-PC ratio # PC
	P01 2 10.1	3.4 0.34 280
	P02 4 18.7	8.7 0.47 825
	P03 5 2153.0	39.8 0.02 1685
	P04 5 50.0	45.1 0.90 1491
	P05 5 113.7	19.5 0.17 1276
	P11 2 86.8	27.4 0.32 513
	P12 4 469.7	50.3 0.11 1821
	P13 5 687.3	176.6 0.26 3454
	P14 5 1198.0	276.7 0.23 4239
	P15 5 1190.9	246.3 0.21 2592

⋆ This work was partially supported by the Ministry of Science and Technology of Spain through MICINN Project DPI2008-02153


Agnieszka Prusińska

Ewa Szczepanik

Alexey Tret'yakov

Singular Lagrange Problem in Calculus of Variations. Elements of p-Regularity Theory

Keywords: singular variational problem, necessary condition of optimality, p-regularity, p-factor operator

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Let us consider the following Lagrange problem:

J 0 (x) = t2 t1
F (t, x(t), x ′ (t))dt → min [START_REF] Brezhneva | Optimality conditions for degenerate extremum problems with equality constraints[END_REF] subject to the subsidiary conditions H(t, x(t), x ′ (t)) = 0, Ax(t 1 ) + Bx(t 2 ) = 0 [START_REF] Brezhneva | Corrigendum: Optimality Conditions for Degenerate Extremum Problems with Equality Constraints[END_REF] where x ∈ C 2 n [t 1 , t 2 ], H(t, x(t), x ′ (t)) = (H 1 (t, x(t), x ′ (t)), . . . , H m (t, x(t), x ′ (t))) T ,

H i : R × R n × R n → R, i = 1, . . . , m, F : R × R n × R n → R, t ∈ [t 1 , t 2 ], A, B - n×n matrices, C l n ([t 1 , t 2 ]
-Banach spaces of n-dimensional l-times continuously differentiable vector functions with usual norms.

Let us introduce a mapping

G(x) = H(•, x(•), x ′ (•)) such that G : X → Y , where X = {x(•) ∈ C 2 n [t 1 , t 2 ] : Ax(t 1 )+Bx(t 2 ) = 0}, Y = C m [t 1 , t 2 ]
. It means that G acts as follows G(x)t = H(t, x(t), x ′ (t)). Then the system of equations ( 2) can be replaced by the following operator equation G(x) = 0 Y (or G(x(•)) = 0 Y ). We assume that all the functions and their derivatives in (1)-( 2) are p + 1-times continuously differentiable with respect to the corresponding variables t, x, x ′ . Under these assumptions: G(x) ∈ C p+1 (X), where by C p+1 (X) we mean a set of p + 1-times continuously differentiable mappings on X.

Let us denote λ(t) = (λ 1 (t), . . . , λ m (t)) T , λ(t

)H = λ 1 (t)H 1 + • • • + λ m (t)H m , λ(t)H x = λ 1 (t)H 1x + • • • + λ m (t)H mx , λ(t)H x ′ = λ 1 (t)H 1x ′ + • • • + λ m (t)H mx ′ .
If Im G ′ (x) = Y, where x(t) is a solution to (1)-( 2), then necessary conditions of Euler-Lagrange

F x + λ(t)H x -d dt (F x ′ + λ(t)H x ′ ) = 0 hold. Here, F x , H x , F x ′ , H x ′
are partial derivatives of the functions F (t, x(t), x ′ (t)) and H(t, x(t), x ′ (t)) with respect to x and x ′ , respectively.

In singular (nonregular or degenerate) case when Im G ′ (x) = Y, we can only guarantee that the following equations

λ 0 F x + λ(t)H x - d dt (λ 0 F x ′ + λ(t)H x ′ ) = 0 (3)
hold, where λ 2 0 + λ(t) 2 = 1, i.e. λ 0 might be equal to 0, and then we have not constructive information of the functional F (t, x(t), x ′ (t)).

Example 1. Consider the problem

J 0 (x) = 2π 0 (x 2 1 (t) + x 2 2 (t) + x 2 3 (t) + x 2 4 (t) + x 2 5 (t))dt → min (4)
subject to

H(t, x(t), x ′ (t)) = = x ′ 1 (t) -x 2 (t) + x 2 3 (t)x 1 (t) + x 2 4 (t)x 2 (t) -x 2 5 (t)(x 1 (t) + x 2 (t)) x ′ 2 (t) + x 1 (t) + x 2 3 (t)x 2 (t) -x 2 4 (t)x 1 (t) -x 2 5 (t)(x 2 (t) -x 1 (t)) = 0, (5)
x i (0)x i (2π) = 0, i = 1, . . . , 5.

Here

F (t, x(t), x ′ (t)) = x 2 1 (t) + x 2 2 (t) + x 2 3 (t) + x 2 4 (t) + x 2 5 (t), A = -B = I 5
, where I 5 is the unit matrix of size 5 and

G(x) = = x ′ 1 (•) -x 2 (•) + x 2 3 (•)x 1 (•) + x 2 4 (•)x 2 (•) -x 2 5 (•)(x 1 (•) + x 2 (•)) x ′ 2 (•) + x 1 (•) + x 2 3 (•)x 2 (•) -x 2 4 (•)x 1 (•) -x 2 5 (•)(x 2 (•) -x 1 (•)) = 0.
The solution of (1)-( 2) is x(t) = 0. At this point G ′ (0) is singular. Later we explain this in more details. The corresponding Euler-Lagrange equation (see [START_REF] Belash | Methods for solving degenerate problems[END_REF]) in this case is as follows:

2λ 0 x 1 + λ 2 -λ ′ 1 + λ 1 x 2 3 + λ 1 x 2 5 -λ 2 x 2 5 -λ 2 x 2 4 = 0 2λ 0 x 2 -λ 1 -λ ′ 2 + λ 1 x 2 4 + λ 2 x 2 3 -λ 1 x 2 5 -λ 2 x 2 5 = 0 2λ 0 x 3 + 2λ 1 x 1 x 3 + 2λ 2 x 2 x 3 = 0 (6) 2λ 0 x 4 + 2λ 1 x 2 x 4 -2λ 2 x 1 x 4 = 0 2λ 0 x 5 -2λ 1 x 5 x 1 -2λ 1 x 2 x 5 -2λ 2 x 2 x 5 + 2λ 2 x 1 x 5 = 0 λ i (0) -λ i (2π) = 0, i = 1, 2.
(to simplify formulas we omit dependence of t here and further in the paper).

If λ 0 = 0 we obtain the series of spurious solutions to the system (4)-( 5):

x 1 = a sin t, x 2 = a cos t, x 3 = x 4 = x 5 = 0, λ 1 = b sin t, λ 2 = b cos t, a, b ∈ R.

Elements of p-regularity theory

Let us recall the p-order necessary and sufficient optimality conditions for degenerate optimization problems (see [START_REF] Brezhneva | Optimality conditions for degenerate extremum problems with equality constraints[END_REF]- [START_REF] Prusińska | P -order Necessary and Sufficient Conditions for Optimality in Singular Calculus of Variations[END_REF]):

min ϕ(x) (7) subject to f (x) = 0, (8)
where f : X → Y and X, Y are Banach spaces, ϕ : X → R, f ∈ C p+1 (X), ϕ ∈ C 2 (X) and at the solution point x of ( 7)-(8) we have:

Im f ′ (x) = Y i.e. f ′ (x) is singular.
Let us recall the basic constructions of p-regularity theory which is used in investigation of singular problems.

Suppose that the space Y is decomposed into a direct sum

Y = Y 1 ⊕ . . . ⊕ Y p , (9)
where

Y 1 = Im f ′ (x), Z 1 = Y.
Let Z 2 be closed complementary subspace to Y 1 (we assume that such closed complement exists), and let P Z2 : Y → Z 2 be the projection operator onto Z 2 along Y 1 . By Y 2 we mean the closed linear span of the image of the quadratic map

P Z2 f (2) (x)[•] 2 .
More generally, define inductively,

Y i = span Im P Zi f (i) (x)[•] i ⊆ Z i , i = 2, . . . , p -1,
where Z i is a chosen closed complementary subspace for (Y 1 ⊕ . . . ⊕ Y i-1 ) with respect to Y, i = 2, . . . , p and

P Zi : Y → Z i is the projection operator onto Z i along (Y 1 ⊕ . . . ⊕ Y i-1
) with respect to Y, i = 2, . . . , p. Finally, Y p = Z p . The order p is chosen as the minimum number for which (9) holds. Let us define the following mappings

f i (x) = P i f (x), f i : X → Y i i = 1, . . . , p,
where

P i := P Yi : Y → Y i is the projection operator onto Y i along (Y 1 ⊕ . . . ⊕ Y i-1 ⊕ Y i+1 ⊕ . . . ⊕ Y p ) with respect to Y, i = 1, . . . , p. Definition 1 The linear operator Ψ p (x, h) ∈ L(X, Y 1 ⊕ . . . ⊕ Y p ), h ∈ X, h = 0 Ψ p (x, h) = f ′ 1 (x) + f ′′ 2 (x)h + . . . + f (p) p (x)[h] p-1 ,
is called the p-factor operator.

Definition 2 We say that the mapping f is p-regular at x along an element h,

if Im Ψ p (x, h) = Y.
Remark 1 The condition of p-regularity of the mapping f (x) at the point

x along h is equivalent to Im f (p) p (x)[h] p-1 • Ker Ψ p-1 (x, h) = Y p , where Ψ p-1 (x, h) = f ′ 1 (x) + f ′′ 2 (x)h + . . . + f (p-1) p-1 (x)[h] p-2 Definition 3
We say that the mapping f is p-regular at x if it is p-regular along any h from the set

H p (x) = p k=1 Ker k f (k) k (x) \ {0}, where Ker k f (k) k (x) = {ξ ∈ X : f (k) k (x)[ξ] k = 0}. is k-kernel of the k-order mapping f (k) k (x)[ξ] k . For a linear surjective operator Ψ p (x, h) : X → Y between Banach spaces we denote by {Ψ p (x, h)} -1 its right inverse. Therefore {Ψ p (x, h)} -1 : Y → 2 X and we have {Ψ p (x, h)} -1 (y) = {x ∈ X : Ψ p (x, h)x = y} . We define the norm of {Ψ p (x, h)} -1 via the formula {Ψ p (x, h)} -1 = sup y =1 inf{ x : x ∈ {Ψ p (x, h)} -1 (y)}. We say that {Ψ p (x, h)} -1 is bounded if {Ψ p (x, h)} -1 < ∞.

Definition 4

The mapping f is called strongly p-regular at the point x if there exists γ > 0 such that sup

h∈Hγ {Ψ p (x, h)} -1 < ∞ where H γ = h ∈ X : f (k) k (x)[h] k Y k ≤ γ, k = 1, . . . , p, h = 1 .

Optimality conditions for p-regular optimization problems

We define p-factor Lagrange function

L p (x, λ, h) = ϕ(x) + p k=1 f (k-1) k (x)[h] k-1 , λ , where λ ∈ Y * , f (0) 1 (x) = f (x) and Lp (x, λ, h) = ϕ(x) + p k=1 2 k(k + 1) f (k-1) k (x)[h] k-1 , λ .
Let us recall the following basic theorems on optimality conditions in nonregular case.

Theorem 1 (Necessary and sufficient conditions for optimality) (see [START_REF] Brezhneva | Optimality conditions for degenerate extremum problems with equality constraints[END_REF]) Let X and Y be Banach spaces, ϕ ∈ C 2 (X), f ∈ C p+1 (X), f : X → Y, ϕ : X → R. Suppose that h ∈ H p (x) and f is p-regular along h at the point x. If x is a local solution to the problem (7)-( 8) then there exist multipliers,

λ(h) ∈ Y * such that L ′ px (x, λ(h), h) = 0 ⇔ ϕ ′ (x) + f ′ 1 (x) + • • • + f (p) p (x)[h] (p-1) * λ(h) = 0. (10)
Moreover, if f is strongly p-regular at x, there exist α > 0 and a multipliers λ(h) such that (10) is fulfilled and Lpxx (x, λ(h), h)[h] 2 ≥ α h 2 for every h ∈ H p (x), then x is a strict local minimizer to the problem (7)-(8).

For our purposes, the following modification of Theorem 1 will be useful (see [START_REF] Belash | Methods for solving degenerate problems[END_REF]).

Theorem 2 Let X and Y be Banach spaces, ϕ ∈ C 2 (X), f ∈ C p+1 (X), f : X → Y, ϕ : X → R, h ∈ H p (x)
, and f is p-regular along h at the point x. If x is a solution to the problem (7)-(8), then there exist multipliers λi (h) ∈ Y * i , i = 1, . . . , p such that

ϕ ′ (x) + (f ′ (x)) * λ1 (h) + . . . + f (p) (x)[h] (p-1) * λp (h) = 0, (11)
and

f (k) (x)[h] (k-1) * λi (h) = 0, k = 1, . . . , i -1, i = 2, . . . , p. (12)
Moreover, if f is strongly p-regular at x, there exist α > 0 and multipliers λi (h), i = 1, . . . , p such that (11)-( 12) hold, and

ϕ ′′ (x) + 1 3 f ′′ (x) λ1 (h) + . . . + 2 p(p + 1) f (p+1) (x)[h] p-1 λp (h) [h] 2 ≥ ≥ α h 2 ,
for every h ∈ H p (x), then x is a strict local minimizer to the problem (7)-(8).

Proof.

We need to prove only the formula (12). From (10) we obtain

ϕ ′ (x) + P 1 f ′ (x) + • • • + P p f (p) (x)[h] (p-1) * λ(h) = 0.
This expression can be transformed as follows

ϕ ′ (x) + f ′ (x) * P * 1 λ(h) + • • • + f (p) (x)[h] (p-1) * P * p λ(h) = 0. Let λi (h) := P * i λ(h), i = 1, . . . , p. Then, for k < i, i = 1, . . . , p, f (k) (x)[h] (k-1) * λi (h) = f (k) (x)[h] (k-1) * P * i λ(h) = = P i f (k) (x)[h] (k-1) * λ(h) = 0, which proves (12).
Now we are ready to apply this theorem to singular calculus of variations problems. Let us introduce p-factor Euler-Lagrange function

S(x) = F (x) + λ(t), g 1 (x) + g ′ 2 (x)[h] + . . . + g (p-1) p (x)[h] p-1 = = F (x) + λ(t)G (p-1) (x)[h] p-1 ,
where

G (p-1) (x)[h] p-1 = g 1 (x) + g ′ 2 (x)[h] + • • • + g (p-1) p (x)[h] p-1
, λ(t) = (λ 1 (t), . . . , λ m (t)) T and g k (x), for k = 1, . . . , p are determined for the mapping G(x) similarly like f k (x), k = 1, . . . , p for the mapping f (x), i.e.

g k (x) = P Y k G(x), k = 1, . . . , p. Denote g (k-1) k (x)[h] k-1 = i+j=k-1 C i k-1 g (k-1) kx i (x ′ ) j (x)h i (h ′ ) j , k = 1, . . . , p,
where g

(k-1)

kx i (x ′ ) j (x) = g (k-1) k x . . . x i x ′ . . . x ′ j (x).
Definition 5 We say that the problem (1)-( 2) is p-regular at x along

h ∈ p k=1 Ker k g (k) k (x), h = 0 if Im g ′ 1 (x) + . . . + g (p) p (x)[h] p-1 = C m [t 1 , t 2 ].
The following theorem holds. k (x). Then there exists a multiplier λ(t) = ( λ1 (t), . . . , λm (t)) T such that the following p-factor Euler-Lagrange equation

S x (x) - d dt S x ′ (x) = F x (x) + + λ, p k=1 i+j=k-1 C i k-1 g (k-1) x i (x ′ ) j (x)h i (h ′ ) j x - ( 13
) - d dt   F x ′ (x) + λ(t), p k=1 i+j=k-1 C i k-1 g (k-1) x i (x ′ ) j (x)h i (h ′ ) j x ′   = 0 holds.
The proof of this theorem is very similar to the one of analogous result for the singular isoperimetric problem, see in [START_REF] Belash | Methods for solving degenerate problems[END_REF], [START_REF] Korneva | Proceedings of Simulation and Analysis in Problems of Decision-Making Theory[END_REF].

Consider again the Example 1 and ( 4)- [START_REF] Prusińska | P -order Necessary and Sufficient Conditions for Optimality in Singular Calculus of Variations[END_REF]. Here p = 2, x = 0. At the beginning we substantiate that G is singular at the points x = (a sin t, a cos t, 0, 0, 0) T .

Indeed, G ′ (x) = (•) ′ 1 -(•) 2 (•) ′ 2 + (•) 1
, where

G ′ (x)x(t) = x ′ 1 (t) -x 2 (t) x ′ 2 (t) + x 1 (t) . Let us denote x ′ 1 -x 2 x ′ 2 + x 1 by x ′ + Lx, where L = 0 -1 0 0 0 1 0 0 0 0 . Then G ′ (x) = (•) ′ + L(•) and
KerG ′ (x) = span {(Φ 1 (t), 0, 0, 0) T , (Φ 2 (t), 0, 0, 0) T } ⊕ {(0, 0, x 3 (t), x 4 (t), x 5 (t)) T ,

x i ∈ C 2 [0, 2π], i = 3, 4, 5}
, where Φ 1 (t) = (sin t, cos t) T , Φ 2 (t) = (cos t, -sin t) T , and moreover Im

G ′ (x) = (Ker(G ′ (x) * ) ⊥ = Ker(-d dt (•) ′ + L T (•)) ⊥ = = ξ ∈ C 2 [0, 2π] : ξ, ψ i = 0, i = 1, 2, ψ 1 (t) = (sin t, cos t) T , ψ 2 (t) = = (cos t, -sin t) T = C 2 [0, 2π].
It means that the mapping G(x) is non-regular at the points x. From the last relation we obtain that

Y 2 = (Im G ′ (x)) ⊥ = span {ψ 1 , ψ 2 } where ψ ′ 1 = ψ 2 , ψ ′ 2 = -ψ 1 and Φ i , ψ j = δ ij , ζ, η = 2π 0 ζ(τ )η(τ )dτ.
The projection operator P Y2 is defined as

P 2 y 1 y 2 = P 2 y = ȳ1 ψ 1 + ȳ2 ψ 2 ,
where y = (y 1 , y 2 ) T and

y -(ȳ 1 ψ 1 + ȳ2 ψ 2 ), ψ 1 = 0, y -(ȳ 1 ψ 1 + ȳ2 ψ 2 ), ψ 2 = 0, i.e. 1 2π y, ψ 1 = ȳ1 , 1 2π y, ψ 2 = ȳ2 . Let us point out that P 2 (x 1 , ψ 1 + x 2 ψ 2 ) = x 1 ψ 1 + x 2 ψ 2 .
Based on Remark 1 we can verify surjectivity of P 2 G ′′ (x)h only on KerG ′ (x), for h ∈ KerG ′ (x) ∩ Ker 2 P 2 G ′′ (x), h = (a sin t, a cos t, 1, 1, 1) T . In order to find P 2 G ′′ (x)h let us determine

G ′′ (x) =           0 0 sin t 0 0 0 0 0 cos t 0 0 0 0 0 cos t -sin t     0 0 cos t 0 0 0 0 0 -sin t 0 0 0 0 0 sin t -cos t           and G ′′ (x)h = 2a 0 0 h 3 sin t h 4 cos t h 5 (cos t -sin t) 0 0 h 3 cos t -h 4 sin t h 5 (sin t -cos t) .
It is obvious that h = (a sin t, a cos t, 1, 1, 1) T belongs to KerG ′ (x) ∩ Ker 2 G ′′ (x) and consequently belongs to KerG ′ (x) ∩ Ker 2 P 2 G ′′ (x). We have

G ′′ (x)[h, x] = 2a(x 3 -x 5 ) sin t cos t + 2a(x 4 -x 5 ) cos t -sin t .
It means that

P 2 G ′′ (x)[h, x] = G ′′ (x)[h, x] and G ′′ (x)[h] • KerG ′ (x) = span {Φ 1 , Φ 2 } = Y 2 .
Therefore G ′′ (x)[h] is surjection. Hence, G(x) is 2-regular along h at the points x = (a sin t, a cos t, 0, 0, 0) T . Finally, we can apply Theorem 3 with λ 0 = 1. We have constructed operator

G ′ (x) + P Y2 G ′′ (x)h = = (•) ′ 1 (•) ′ 2 +
0 -1 2a sin t 2a cos t 2a(cos t -sin t) 1 0 2a cos t -2a sin t 2a(sin t -cos t)

which corresponds to the following system (F x ′ = 0):

F x (x) + (G ′ (x) + P 2 G ′′ (x)h) * λ = 0 ⇔ F x (x) + G ′ (x) T λ + (P 2 G ′′ (x)h) T λ = 0, ⇔ ⇔            2x 1 -λ ′ 1 + λ 2 = 0 2x 2 -λ ′ 2 + λ 1 = 0 2x 3 + 2λ 1 a sin t + 2λ 2 a cos t = 0 2x 4 + 2λ 1 a cos t -2λ 2 a sin t = 0 2x 5 + 2λ 1 a(cos t -sin t) + 2λ 2 a(sin t -cos t) = 0 (14) or ⇔                2a sin t -λ ′ 1 + λ 2 = 0 2a cos t -λ ′ 2 + λ 1 = 0 λ 1 sin t + λ 2 cos t = 0 λ 1 cos t -λ 2 sin t = 0 λ 1 (cos t -sin t) + λ 2 (sin t -cos t) = 0, λ i (0) -λ i (2π) = 0, i = 1, 2.
One can verify that the false solutions of (6)

x 1 = a sin t, x 2 = a cos t, x 3 = x 4 = x 5 = 0 do not satisfy the system (14) for a = 0. It means that x 1 = a sin t, x 2 = a cos t, x 3 = x 4 = x 5 do not satisfy 2-factor Euler-Lagrange equation (13) Let us consider the same problem with higher derivatives x ′ (t), . . ., x (r) , r ≥ 2,

J(x) = t2 t1 F (t, x(t), x ′ (t), . . . , x (r) (t))dt → min, x(t) ∈ C 2r n [t 1 , t 2 ],
subject to subsidiary differential relation

H(t, x(t), x ′ (t), . . . , x (r) (t)) =   H 1 (t, x(t), x ′ (t), . . . , x (r) (t)) • • • H m (t, x(t), x ′ (t), . . . , x (r) (t))   =   0 . . . 0   , A k x (k) (t 1 ) + B k x (k) (t 2 ) = 0, where A k , B k are n × n matrices, k = 1, . . . , r. Let G(x) = H(•, x(•), . . . , x (r) (•)), G : X → Y, where Y = C m ([t 1 , t 2 ]) and X = {x(•) ∈ C 2r n [t 1 , t 2 ] : A k x (k) (t 1 ) + B k x (k) (t 2 ) = 0, k = 1, . . . , r}.
Moreover,

g (k-1) k (x)[h] k-1 = i1+•••+ir=k-1 g (k-1) x i 1 •••(x (r) ) ir [h + h ′ + • • • + h (r) ] k-1 , k = 1, . . . , p,
and introduce the co called p-factor Euler-Poisson function

K(x) = F (x) + λ(t), g 1 (x) + g ′ 2 (x)[h] + . . . + g (p-1) p (x)[h] p-1
Theorem 4 Let x(t) be a solution of the problem (1)-( 2) and assume that this problem is p-regular at x along h ∈ p k=1 Ker k g (k)

k (x). Then there exist a multiplier λ(t) = ( λ1 (t), . . . , λm (t)) T such that the following p-factor Euler-Lagrange equation

K x (x) - d dt K x ′ (x + d 2 dt 2 K x ′′ (x) -... + (-1) r K x (r) (x) = = F x (x) + λ(t), p k=1 g (k-1) k (x)[h] k-1 x - - d dt F x ′ (x) + λ(t), p k=1 g (k-1) k (x)[h] k-1 x ′ + + . . . + (-1) (r) d r dt r F x (r) (x) + λ(t), p k=1 g (k-1) k (x)[h] k-1 x (r) = 0
holds.

The proof of Theorem 4 is similar to that one the reader can find in [START_REF] Korneva | Proceedings of Simulation and Analysis in Problems of Decision-Making Theory[END_REF] for isoperimetric problem. Example 2. Consider the following problem

J(x) = π 0 (x 2 1 (t) + x 2 2 (t) + x 2 3 (t))dt → min (15)
subject to

H(t, x(t), x ′ (t), x ′′ (t)) = x ′′ 1 (t) + x 1 (t) + x 2 2 (t)x 1 (t) -x 2 3 (t)x 1 (t) = 0, (16)
x i (0) -x i (π) = 0, x ′ i (0) + x ′ i (π) = 0, i = 1, 2, 3. Here A 1 = -B 1 = I 3 , A 2 = B 2 = I 3
, where I 3 means the unit matrix of size 3.

The solution of (15)-( 16) is x(t) = 0. The Euler-Poisson equation in this case has the following form

λ 0 F x + λ(t)H x - d dt (λ(t)H x ′ ) + d 2
dt 2 (λ(t)H x ′′ ) = 0 or 2λ 0 x 1 + λ + λx 2 2λx 2 3 + λ ′′ = 0 2λ 0 x 2 + 2λx 2 x 1 = 0 2λ 0 x 3 -2λx 3 x 1 = 0, λ(0)λ(π) = 0, λ ′ (0) + λ ′ (π) = 0 and gives us the series of spurious solutions x 1 = a sin t, x 2 = 0, x 3 = 0, λ = b sin t, λ 0 = 0, a ∈ R. The mapping G(x) is singular at these points x 1 = a sin t, x 2 = 0, x 3 = 0 and G ′ (a sin t, 0, 0) is non surjective. But G(x) is 2-regular at the points x = (a sin t, 0, 0) along h = (sin t, sin t, -sin t). Indeed, Y 2 = span {sin t}, G ′ (x)h + P Y2 G ′′ (x)[h] 2 = h ′′ + h + 2a sin t The first equation has no solutions for a = 0, which means that the point x = (a sin t, 0, 0) T is not a local solution of the considered problem.

Theorem 3

 3 Let x(t) be a solution of the problem (7)-(8) and assume that the problem is p-regular at x along h ∈

π 0 ( 0 (

 00 sin 2 t -sin 2 t) sin 2 tdt = 0.It means that h ∈ KerG ′ (x) ∩ P Y 2 Ker 2 G ′′ (x) andP Y2 G ′′ (x)h = 2a sin t π sin t(•) 2 + sin t(•) 3 ) sin 2 tdt.We haveP Y2 G ′′ (x)h   b sin t b sin t b sin t   = 2ab sin t π 0 2 sin 4 tdt = Y 2 , b ∈ R i.e.G is 2-regular at the points x = (a sin t, 0, 0) along h. At these points x we can guarantee λ 0 = 1 in the 2-factor Euler-Poisson equation 2a sin t + λ ′′ + λ = 0 2a sin t π 0 sin 3 τ λ(τ )dτ = 0 2a sin t π 0 sin 3 τ λ(τ )dτ = 0 λ(0)λ(π) = 0, λ(0) + λ(π) = 0

Acknowledgements. Research of the second author is supported by the Russian Foundation for Basic Research Grant No 11-01-00786a and the Council for the State Support of Leading Scientific Schools Grant 5264.2012.1.