Dr Gérald Billard

Ludivine Lemoine

THE GEOGRAPHY OF ICE HOCKEY IN FRANCE : THE WEIRD LOCATION OF PROFESSIONAL TEAMS

The lock-out affecting the NHL since September 16 th 2004 had two consequences. Firstly, five months after the lock-out began, the NHL commisioner Garry Bettman cancelled the season in February 2005. Secondly, around 180 NHL players have left North-America for European teams. Thus 2 well-known players from the Calgary Flames (NHL 2004 Stanley Cup finalists) Steve Montador and Steven Reinprecht have joined the team of Mulhouse. This arrival, and others like in Villard de Lans (French Alps), had enlightened professional ice hockey in France. It was a great opportunity for this minor sport in France to be under the light of media and gain more consideration.

According to the French Minister of Sports, 17 229 people were practicing ice hockey in France in 2004 out of 60 millions inhabitants. This is far from the two millions French people who play soccer 1 , one million tennis, 564 000 judo or 426 888 basket-ball. Professional soccer is also the most attended sport on television in France. In 2004, the ten best audiences, regarding all programs (sports, movies, news magazines…) were made owing to 10 soccer games 2 . And the (private) TV channel Canal+ broke a new record in 2004, by making an agreement with the Professional Soccer League 3 for 600 millions euros per year (for three years) to have the right to broadcast the French first league championship. A contrario, the professional hockey games have just started to be on air 4 : this year, a cable and satellite TV channel Sport +, is broadcasting only few games of the regular season, the Magnus Cup 5 finals series and the French All Stars Game.

But, compared to other ice sports, like figure or speed ice skating, ice hockey has more permit-holders. After few years of confusion at the end of the 90's, the professional hockey league seems to have found a sustainable formula and seduced the fans. In some cities, like Rouen, Amiens and Grenoble the mean attendance per game rises over 2 000 spectators 6 . It is interesting to point out that professional ice hockey has not been much investigated by the French geographers 7 despite the captivate problematic of the spatial diffusion of professional ice hockey in France.

1 20 % of the total number of players in all sports in France play soccer and only 0.38 % play ice hockey 2 Eurodata TV Worldwide, march 2005 3 Ligue Professionnelle de Football (LFP) 4 In 2003, the Magnus Cup finals series were not broadcasted on television , only games 2 and 3 of the finals series were broadcasted live on an Internet site (www.planisport.com). 5 Magnus Cup is the French name of the professional trophy for the ice hockey, like Stanley Cup for the NHL. 6 Source : www.hockeyarchives.

The slow diffusion of ice hockey in France

Most of the leading sports in France have their global federation and a specific association for professional teams, like the French Federation of Soccer and the Professional Soccer League. Ice hockey does not have a dedicated association and shares a federation with the other ice sports (Skating, Bobsleigh..). It is obvious that ice hockey is a minor sport in France despite the fact that the French ice hockey championship celebrated its centenary in 2003.

At the beginning, ice hockey was a seasonal sport played only in winter and outside. Paradoxically, in France, ice hockey did not appear in the mountains (Alps, Pyrénées…) or in the coldest regions but in Paris, on an iced canal in Versailles by the year 1891. Ten years later, Pierre de Coubertin, the father of the modern olympics games, created the first French ice hockey team : the Paris Hockey Club. They won against the London Princess Ice Hockey team which was seven times world champion. In 1905, several teams were created in Paris (and one in Lyon, first champion that same year) and they dominated the French league until 1922 (Biot, 1996).

In the 1920's, a new pole of spreading appeared in the French Alps. Many teams were created in several cities (Chamonix, Gap..) without any logic regarding the size of the urban areas. Grenoble, the biggest town in the Alps, had its own team later than several smaller mountain towns. All the cities in the Alps share an ideal climatic environment to host hockey players with a cold and long winter which is essential to practice ice hockey outdoor during a long period. Chamonix became the major French team ever by winning 30 championship titles between 1923 and 1979 (see Fig. 1).

In 1930, Chamonix hosted for the first time the Ice Hockey World Championship. Nevertheless, it is interesting to note that ice hockey was not present in the other French moutains (Massif Central, Pyrénées, Jura…).

After the 1968 Olympics Winter Games in Grenoble, ice hockey became more popular in France and spreaded in the northern part of the country, in the plains. Some cities such as Amiens, Rouen, Viry-Châtillon, Angers or Reims built their first indoor ice rink. In the 70's and the 80's, with the construction of a lot of indoor ice rinks in the large cities, several new teams were created. This spreading of ice hockey was enabled by the large investments made by the local governments, especially the local authorities.

The modern structuration of French professional ice hockey

With the television covering of the French team games8 during the 1992 Olympics Winter Games in Alberville (French Alps), more and more young people decided to play ice hockey. In France, between 1992 and 1993, the number of players rised from 29 000 to 32 000. But, the best professional teams were affected by a major financial crisis started in 1991. Only six teams played the 1991-1992 championship (Magnus Cup) and ten in the minor division.

In 1992-1993, a lot of teams (Briançon, Bordeaux, Epinal or Clermont-Ferrand) had not enough money to compete in the Magnus Cup again. Errors in the financial management of the teams, the retreat of big private sponsors (ex. Bayer Inc. in Caen), less money from the local authorities could explain this crisis. Only four teams, survivors from this major failure, joined the new Magnus Cup formula composed of 16 teams.

We can also argue that despite the good results of our national team, there is a lack of a real ice hockey culture in our country illustred by the weak coverage of ice hockey events by the media. Furthermore, ice hockey does not run its own federation and has to share the French Ice Sport Federation9 (FFSG) with some traditional non-professional sports like bobsleigh or curling, ! For this reason, it is very difficult to develop professional skill.

The season of 2003/2004 seemed to be very crucial for the professional ice hockey because it negociated (and gained) more self-sufficiency from FFSG and signed some agreement with a national television channel for the broadcasting of the Magnus Cup Finals. Since 2004, the sitution of the teams has improved owing to the creation of the French Hockey Executive Authority (AEHF), under the direction of a former pro player (Luc Tardif). This authority wanted to bring some cohesion back in ice hockey as well as a better management of the professional teams. In 2003, a first rule was settled to restrict the total salaries of each professional team : 450 000 €/year was the new limit. Then in june 2004, only 4 non-european players could be hired in professional teams of the Super 16 (see Fig. 2). All these decisions show the will of the AEHF to have a sane financial management in professional ice hockey and encourage the training of French young players instead of importing, as usual, canadians ones, for example.

The weird location of the French champion teams in the history of ice hockey

If the formula of the Magnus Cup has changed many times over the last century, the results (and sometimes banckruptcy !) always determine if a team will stay in Magnus Cup division or will be replaced by a team of the lower level division. In the NHL, the system of professional franchises frozes the location of professional teams (except sale and relocation of franchises). In France, the geography of the distribution of professional teams should be uncertain in times, even it would be logical to find more good teams in the coldest and icy regions like in the mountain areas.

Nevertheless, the first concentration of ice hockey teams was in Paris area (17 Magnus Cup titles won between 1903 and 1939) and in Lyon, the two largest cities in France, both located in the plain. In the begining of the twentieth century, this phenomena could be explained by the role traditionally played by large urban areas in the diffusion of innovation. Until 1922, Chamonix (see Fig. 3) was the first town in the Alps regions to have won a Magnus Cup : that reaches our first hypothesis that weather conditions could be a good factor to understand the development of ice hockey (Biot, 1996). In fact, since 1990, 15 Magnus Cups have been attributed and only Grenoble (Alps region) won the ultimate title in 1991 and 1998 over the big domination of the cities located in the plains (Rouen, Amiens, Reims) or on the coastline (Brest) 10 . The teams in Paris and Lyon disappeared from the Super 16 division. More globaly, the current analysis of the development of the professional teams in France shows a triple phenomena : an historical concentration in the Alps region, a large diffusion of outstanding teams in the plain areas around Paris in the western and eastern parts of France (north of the Loire river) and a lack of professional teams in the centre region (Massif Central) and in southern regions (French Riviera, South-West), except in the city of Anglet. (see Fig. 2).

The primary factors of location of the professional teams in France

Since the end of the 60's, cities located in the northern part of France, in the plain area, have been experimenting the creation of several professional teams. This diffusion has had an impact on the geographical distribution of teams able to win the Magnus Cup and we wonder which factors could explain this new spatial concentration of the best French professional teams. The location of the professional teams shows a country cut in two parts along a line between Paris and Lyon. At the beginning of the XX th century, players were practising only outdoor and in winter. If Paris is the first place in the development of ice hockey, we can point out the influent role played by Pierre de Coubertin. Then, the domination of the team of Chamonix on French ice hockey can be explained by weather conditions : with long and cold winters, the outdoor ice rinks could be built for several months. Moutains are also a good place for the development of a winter sport culture.

Since 1955, the massive construction of indoor ice rinks in France (see Fig. 4), especially in the large cities located in the northern plains, has changed the geography of ice skating and since the 60's, of the creation of ice hockey teams. We can link the construction of indoor sport facilities with the climatic factor because indoor sports (hand-ball, volley-ball, basket-ball…) are more common on the northern regions of France whereas outdoor sports (tennis, rugby, surf…) are more frequent in the South. Also, we can argue that cities in the South are far away from the historical first poles of ice hockey (Paris and the Alps) and maybe, it explains the weakness of the ice hockey culture in the southern regions too (Lemoine, 2001).

Another element can help us to understand the geography of ice hockey in France : the global concentration of sport facilities construction in the most dynamic regions in terms of population and economic growth. Indeed, the construction of an ice hockey stadium is a big financial weight for the local communities. Lemoine (2001) shows a correlation between high density urban areas and an important number of ice hockey professional teams. For example, since 1960, the western part of France has been experimenting a large growth of population and hockey teams as well. On the contrary, the Massif Central or the south-west region have a small density of population and only few ice hockey teams. But the number of inhabitants is not the only factor to be considered. An ice hockey rink is not an usual facility for a city on the contrary to gymnasiums, swimming pools, soccer or rugby stadiums or tennis courts. For example, in 1998, there were 3 000 indoor and 10 000 outdoor ice rinks in Canada (for 30 millions inhabitants), 1 500 and 2 500 in the United States but only 97 indoor and 48 outdoor ice rinks in France (Evdokinoff, 1998). Ice hockey is one of the only sports in France that can not be played outdoor because of the temperated weather and anywhere (outside an ice rink). The financial cost of building and maintaining an ice rink is so high that only large cities or winter stations with a lot of tourists can afford such an unusual sport facility. Furthermore, the level of human density is an important factor to reach the financial rentability of an ice rink, to generate enough skaters and potential hockey players.

The case of Rouen's Dragon, to seize the complexity of the development of ice hockey in France

In a previous part of this paper, several structural elements have been identified as primary factors to explain the spread of professional ice hockey in France. Nevertheless, to understand the complexity of the birth and the rise of a professional team at the ultimate level, we have to go from the macro to the micro scale. Thus, the study of the development of ice hockey in the city of Rouen (400 000 inhabitants in the metropolitan area) located in Normandy is very interesting to point out the numerous factors with which we can try to understand a success winter sport story in a western temperate region.

The first days of ice hockey in Rouen

Like in many cases, the development of a new sport in a determinated place is often a question of luck and goodwill of one or several human beings. In 1967, Jean-Jacques Calvé, decided to build a temporary ice rink on the riverside of the Seine for leisure skating.

Then, Bernard Le Feuvre, who had worked in Zurich (Switzerland) with a former hockey player and François Legay, found of ice hockey after several stays in the moutain areas, introduced ice hockey to the people of Rouen for the first time (Biot, 1996 ). Many young inhabitants were interested in this sport but a lack of basic equipments (stick, gloves, helmets, skates...), the few hours dedicated to ice hockey practising and the competition with the soccer were negative elements that frozed ice hockey expansion. But four years later, the City of Rouen decided to build a modern ice rink with a capacity of 876 seats. When the ice rink was opened in 1971, the Rouen Olympic Club (ROC) was already born (1970), composed of both skaters and hockey players. At first, hockey players had lots of problems for practising because they could not use the ice more than one hour and fifteen minutes twice a week and they did not have the right to play with a puck because it made some damages to the white paint around the ice rink !

In November 1972, the first hockey team of Rouen played its first game and lost 22 to 0 against Caen. The first official season for the ROC was in 1973/1974 and two foreign players joined the team. Many inhabitants were curious about this very new physical and intensive sport. In 1978, the birth of the Rouen Hockey Club (RHE) was an important step because the hockey players had for the first time their own club (without ice figure skaters). The 80's were very important for the RHE which had 65 players with two hours of practice a week in 1979. Then it became an important club of 350 players with 7 teams in different championships and a school of hockey for kids. In competition with ice figure skating, hockey had not enough time for praticing on the ice but, according to Lemoine ( 2001), local authorities gave 1 850 € each year to the hockey team (only 540 € in the 70's).

The arrival of a new trainer and new players, especially from Canada or northern countries of Europe (Sweeden, Finland…), and the setting on of a new politic in the training of the young players had transformed ice hockey in Rouen. It became the most popular sport because, at the same time, the local soccer professional team declined. Twelve years after the opening of the ice rink, the Rouen Hockey Club got to the professional top national division.

A fan was so happy with the good results that he gave the nickname Dragons11 to the players of Rouen and all the fans, medias and players themselves aggreed : the Rouen Dragons were born. Ice hockey became more and more popular and the mean attendance for each game rised from 400 in 1985 to 1 500 fans in 1986. The need for a permanent financial and organizational management was another step in the way of professionalism.

The factors of the rise of the Dragons

The 90's were years of consecration for the dragons. The Rouen Dragons became the greatest professional team in France with seven Magnus Cup victories from 1990 to 2003 and an European Title in 1996.

In 1988-1989, the arrival of new sponsors increased the economical power of the RHC. A lot of good foreign players arrived and the Dragons went to training camps in Europe. The economic factor is a crucial element for the development of professional hockey.

In France, without television rights, low income from merchandising, and small attendances, the French ice hockey depends on the private sector (local compagnies). The will of RHC to stay on the top of the professional hierarchy (Magnus Cup in 1990, 1992, 1993, 1994, 1995) goes with an invariably inflationist politic : in 1997, RHC was bankrupt because of a too heavy financial pressure It is also an interesting program to build localy a real ice hockey culture.

To conclude, the success of the Dragons can be explained by a complex interactions of several positive and complementary factors (see Fig. 7) which start by the involvment of a community of enthusiastic people (fans, parents of the young players, volunteers…), private sponsors, local governments, set up of a performant and sustainable training school.

In Rouen, the lack of other professional teams in more coverage sports (soccer, basket-ball…) is also an essential factor. In this condition, to understand the development and the location of the best professional ice hockey teams in France, we often conduct our analyses at the local scale and dive in the history of each team.

BALANCE SHEET OF ONE CENTURY OF ICE HOCKEY IN FRANCE

Far away from a large and total territory coverage, ice hockey in France seems to be a minor sport with only two historical poles (Paris and the Alps) and a current development in large cities of the northern part of the country explained by a local complex interaction system. On the contrary to France, ice hockey in Canada is the national sport and is part of the construction of the Canadian Nation and of the construction of a specific identity as regard to the Unites States and Europe (Augustin, 2000).

In 1998, there were 29 679 teams and 465 000 players in Canada and 22 234 teams and 294 195 players of ice hockey in the United States (Evdokinoff, 1998). Very spectacular but paradoxaly with a small media coverage in France, professional ice hockey is subject to a lot of constraints : financial, cultural, climatic, political... Nevertheless, there are 475 teams in France and 10 564 players. But in the 30 professional teams in France, a lot of players must have a second job because of the small salary. Despite a lack of information, we can say that the mean salary for a professional player is around 1 500 €/month : for the rookies, this salary is around 1 000 €/month. The best salaries are often reserved for foreign players. In 2004, Merrigan and Trudel (2004)

Fig. 1

 1 Fig.1 French teams winners of the Magnus Cup between 1903 and 2005

  But beyound the financial aspect, the long-term stability of the Magnus Cup is for us an important goal to gain more professionalism in ice hockey and to develop loyalty of the fans. If the Magnus Cup has been existing since 1903, the championship has changed its formula several times. If we consider the last five years, five formulas of the Magnus Cup have been tested with seven teams during the season 2000/2001 and 15 teams in 2004/2005 ! Now, in the Super 16 (Magnus Cup), 15 teams play in the regular season in a unique division with a first and second leg. The 8 best teams run in play-off (3 wins to be qualified) and at the end, one of the finalists need two victories to win the Magnus Cup. The teams not qualified for the play-offs run a relegation division and the teams with the more losts have to leave the Super 16 to play at a lower level the following year.

Fig. 2

 2 Fig.2 Teams in Super 16 (Magnus Cup) for season 2004/2005

Fig. 3 :

 3 Fig. 3 : Players practising ice hockey outdoor in 1920's in Chamonix (Alps area -France) (Source : Biot R, 1996)

Fig. 4

 4 Fig. 4 The construction of ice rinks in France after 1955

  . The creation of a new team R.H.E. 76 (Rouen Hockey Elite 76), a private compagny with a sport object (S.A.O.S., law n° 84.61 of july 16th 1984) under the control of the AEHF more and more vigilant about the financial management of professional teams 12 : in 2001 and 2003, the Dragons won two more Magnus Cups (see Fig. 5).

Fig. 5

 5 Fig. 5 The Rouen Dragons team of season 2004/2005 (source : www. rouenhockeyelite76.com)

Fig. 6

 6 Fig.6The olympic ice rink located on Lacroix Island, home of the Rouen Dragons (G. Billard ©)

Fig. 7

 7 Fig.7 The complex system of the development of professional ice hockey in the city of Rouen

  pointed out that the mean salary for a NHL player was around 1,7 million / year. Now, out of the 304 players playing in Super 16, 104 are from foreign countries (34 % of the total). This foreign légion empowers French teams which have some difficulties to form good French players. In fact, only two French players have played in the prestigious National Hockey League ever : Philippe Bozon (St Louis Pinguins) and Christobal Huet (Montreal Canadians). La culture du sport au Québec, Centre d'études canadiennes, Talence : MSHA •Augustin, (2000), "La culture sportive au Canada : le hockey sur glace comme sport identitaire", in Sport et identités, Paris, l'Harmattan : 69-72 • Biot R, (1996), Rouen, les années Dragons, Rouen : Sports de Glace Editions • Evdokinoff S., (1998), Le hockey sur glace, Toulouse : éditions Milan • Grelon B., (1998), Cours de hockey sur glace, Paris : éditions De Vecchi • Lemoine L., (2001), Le sport, une culture : l'exemple du hockey sur glace rouennais, Master's thesis, sous la direction de Françoise Lucchini, Université de Rouen. • Le journal du hockey, n°49, avril 2001, Grenoble. • Merrigan P. Trudel P., (2004), Dernières minutes de jeu : les millions du hockey, Paris, Editions Hurtubise HMH.

RÉFÉRENCES

• Augustin J-P, (1995), Sport, géographie et aménagement, Paris : Nathan •Augustin J-P, Sorbets C.,

(1996)

,

After some fantastic games, the French National Team finished at the 8 th position of the Olympic tournament.

Fédération Française des Sports de Glace(FFSG)

Source : www.ffsg.org

The dragon is a mystical animal on the front of the drakkars, boats of the Vikings who invaded Normandy before the Middle Age. A dragon was also the emblem of a valiant knight who battled fidely closed on Guillaume Le Conquérant's side.

In March 2005, the FFSG announced that 4 teams of the Super 16 would not be authorized to join the Magnus Cup 2005/2006 if their financial statement did not improve. (source : www.ffsg.org).

For comparison, the new 375 millions $ Staples Center in Los Angeles, home of the LA Kings, has 19 000 seats for hockey.

40 young players (from 11 to 18 years old) from 2 high schools are involved in this program


Kshitij Kulshreshtha

Sri Hari Krishna Narayanan

Tim Albring

A mixed approach to adjoint computation with algorithmic differentiation

Keywords:

Various algorithmic differentiation tools have been developed and applied to large-scale simulation software for physical phenomena. Until now, two strictly disconnected approaches have been used to implement algorithmic differentiation (AD), namely, source transformation and operator overloading. This separation was motivated by different features of the programming languages such as Fortran and C++. In this work we have for the first time combined the two approaches to implement AD for C++ codes. Source transformation is used for core routines that are repetitive, where the transformed source can be optimized much better by modern compilers, and operator overloading is used to interconnect at the upper level, where source transformation is not possible because of complex language constructs of C++. We have also devised a method to apply the mixed approach in the same application semi-automatically. We demonstrate the benefit of this approach using some real-world applications.

Introduction

Solution techniques for optimal control and optimal design problems rely on the correct and efficient computation of the adjoint state. Various analytical and numerical techniques have been devised to compute these derivatives in the past. One of the emerging techniques for the computation of derivatives on modern computers is algorithmic differentiation (AD) [?]. Despite differentiation being a badly conditioned operation in general, research has shown [?] that the process of algorithmic differentiation is well behaved and the derivatives obtained are accurate to within round-off errors. This situation is in contrast to numerical derivatives computed by using finite-differencing techniques, where the difference step size is of critical importance.

Algorithmic differentiation assumes that functions are evaluated by using a finitely terminating evaluation procedure consisting of simple arithmetic operations {+, -, /, * } and elementary function evaluations { √ , sin, cos, exp, log, . . .}. Since the analytical derivatives of such arithmetic operations and elementary functions are well known, these can be introduced in the evaluation procedure almost mechanically, and the chain rule of differentiation can be applied to propagate the derivatives from one variable to another in the evaluation. In [?] various modes of propagation of derivatives as well as methods to implement tools are discussed in great detail. Here we present two techniques of AD, namely, source transformation and operator overloading.

Source transformation: Source transformation AD tools generate a new source code that computes the derivatives of an input source code. The output code must be compiled and executed in order to compute the derivatives. Tools such as ADIFOR [?], Tapenade [?], and OpenAD [?] can be used to generate derivative code for functions written in Fortran. Tapenade and ADIC [?] are examples of source transformation AD tools for C. In this work, we use ADIC to differentiate input source code portions written in C. In the output code, active variables are declared as objects of DERIV TYPE, and runtime functions are used to propagate derivatives between them. When the output code is compiled with an appropriate driver and runtime library provided by ADIC, the Jacobian matrix can be computed.

Because such tools perform source code analysis, they can identify algorithmically active and passive variables and portions of the code. Furthermore, compilers can optimize the output code, resulting in high performance. However, no tool can generate derivative code for complete C++ input. C++ contains features such as polymorphism, inheritance, and templates that cannot be resolved statically, precluding the generation of correct derivative code.

Operator overloading: In an object-oriented language such as C++ the concept of operator overloading is well known. Several tools have been developed in recent years for AD using C++ operator overloading. ADOL-C [?] is a well-known open source AD tool with many features and high flexibility and has been successfully used to compute derivatives in a large number of simulation codes. The most important manual change required for using ADOL-C in any simulation is to convert the datatype of the real values to the special datatype adouble defined in the ADOL-C library. All operations executed after a call to trace_on() and before a call to trace_off() are recorded in an internal representation called the trace. Before the actual computation takes place, the independent variables are marked by assigning them values using the special <<= operator. Similarly the final dependent variables are marked by extracting their values using the special >>= operator. The trace can then be used in any mode of AD (i.e., forward or reverse) in order to compute first or higher derivatives. Several easy-to-use drivers for computing the derivative information from the trace are available. The most-used drivers are gradient(), jacobian(), and hessian(). For further usage details see [?].

The creation of the trace is the most crucial part of the whole program; and depending on the complexity of the functions being traced, the trace can become large and thus has the most impact on the memory consumption of the program. Where the trace does not fit into a prescribed amount of memory (RAM), it spills over automatically to the disk as trace files, thereby reducing the performance of the implementation severely. Past attempts at reducing the memory requirement for certain applications include using checkpointing strategies [?,?]. For a number of problems, however, checkpointing is not applicable.

We propose a mixed approach that uses both operator overloading and source transformation to differentiate an input code that is largely written in C++ but whose computationally intensive portions are written in a C-like manner. Our approach employs operator overloading for most of the application and source transformation for the C-like portions. Because the computationally intensive portions contribute most to the trace, using source transformation instead for these portions leads to a smaller trace and better performance. We have made changes to both ADIC and ADOL-C and written a preprocessor that enables the approach to be semi-automated. The rest of the paper is organized as follows. Section 2 presents the details of the mixed approach. Section 3 presents experimental results on two applications and, Section 4 discusses future work.

Mixed approach

The process of converting an ADOL-C instrumented application to use ADIC in certain parts is the following: (1) The user identifies a computationally intensive and C-like function from the input based on performance analysis or experience.

(2) This function must be treated as an externally differentiated function (EDF) by ADOL-C. For this purpose, annotations are added to the input to support extraction of the EDF and its callees for differentiation by ADIC. Additional annotations are used to generate wrappers functions and files to copy data between ADOL-C data structures and the EDF. (3) ADIC is used to differentiate the EDF and provide forward-and reverse-mode differentiated code for it. (4) The EDF input, output, wrapper files, and original ADOL-C code are then built together. The rest of this section elaborates on the concepts of the EDF and the changes we made to ADOL-C and ADIC to support the mixed approach.

Externally differentiated functions in ADOL-C:

The individual arithmetic operations and mathematical function evaluations of an EDF are not recorded on the ADOL-C trace. Instead the actual implementation of the differentiated EDF is provided via user defined function pointers that implement a certain predefined signature. As one can see in Fig. 1 the EDF replaces a large part of the trace by repeated calls to itself, which reduces the size of the trace. When ADOL-C processes the trace and arrives at a call to the EDF, ADOL-C calls the corresponding user-provided forward mode-or reverse mode-derivative code to obtain the derivatives.

ADOL-C previously maintained the EDF interface using a special structure struct ext_diff_fct that is registered to the ADOL-C core on a per function basis. Implementations for the forward-and reverse-mode first-order derivative computations are set up in this structure as function pointers that have a particular signature as defined in the header file <adolc/externfcts.h>. The limitation of this interface is that it expects all inputs as well as all outputs to the EDF to be passed as two contiguously allocated arrays. The design of the adouble type in ADOL-C creates an internal representation of the executed code at runtime. In order to do so most efficiently, adouble objects are allocated in a memory pool whereever there is unused space. Unless the pool is exhausted, new memory is not allocated. This design makes the allocation of large contiguous arrays an expensive operation, because of the need for finding a suitable chunk of unused space in the memory pool. Several smaller contiguous arrays, on the other hand, can be allocated more easily. Therefore we designed a second version of the EDF interface structure struct ext_diff_fct_v2 that supports providing several input arrays and several output arrays, each not necessarily of the same size. We also added extra integer-valued input parameters and an opaque object-valued input/output parameter that do not have an effect on the differentiation process outside the EDF. These changes required adjusting the signatures of the forward-and reverse-mode implementations for the EDF. The signatures now contain the number of input and output vectors, the sizes of each of these, their values, the corresponding tangents or adjoints, extra integervalued input arrays, and an opaque context object if needed (see Fig. 2(a)). However, the process of registration and setup of the function pointers stays the same as in the original EDF interface and can even be encapsulated in a separate routine (see Fig. 2(b)), which is called once before the function is required to be evaluated. The ADOL-C evaluation of the complete structure would then look something like the code in Fig. 2(c). The user-provided functions edf-> fov_forward() or edf->fov_reverse() are called during the evaluation of the jacobian() at the appropriate point.

Runtime support for ADIC generated code:

To support the mixed approach's use of forward-and-reverse mode AD in a single execution instance, we recoded ADIC's runtime library in C++ and used namespaces to separate forwardand reverse-mode derivative manipulation routines. The namespace usage is inserted into the ADIC-generated code by using simple postprocessing scripts. The DERIV_TYPE structure was rewritten to be a class that supports both dynamic and static allocation of the grad array within DERIV_TYPE. Because dynamic allocation for every DERIV_TYPE object can be expensive, we created a memory manager that allocates a large amount of memory from the heap and then allocates the User annotations and preprocessing: User annotations have two purposes: First, they identify an EDF and its callees for extraction and subsequent differentiation by ADIC. The annotations surround the EDF and its callees, as shown in Fig. 3(a). The extraction of code is necessary because ADIC requires the C code to be isolated from the C++ code that it does not differentiate. Additional user editing may be required to obtain code that is appropriate for differentiation by ADIC. Second, annotations are used to generate the interface code for arguments of the EDF. The annotation identifies inputs, outputs, and their respective sizes or extra integers required for the computation, as well as the position of each formal parameter in the argument list of the EDF. This information helps generate wrapper code to transfer data between ADOL-C data structures and ADIC-generated code. These annotations are written directly as Python tuples, as seen in Fig. 3(b). Each tuple contains the name of the formal argument, followed by its size and the position in the formal argument list. The size itself is a list of length 0, 1, or 2, depending on whether the argument represents a scalar, a vector, or a matrix. Integer arguments are always scalars. The size may also contain references to values stored in the integers list. Several interface definitions, and thus multiple EDF structures, may also be used in any application.

/*@ declare doubletype=adouble @*/ // since ADIC doesn't know adouble /*@ begin adic_extract global @*/ ... // global defines, variables etc. // required by ADIC routines /*@ end adic_extract @*/ /*@ begin adic_extract @*/ ... // lower level computational routines // differentiated by ADIC /*@ end adic_extract @*/ /*@ begin adic_extract replace=rk_iter type=void @*/ ... // top level interface routine // differentiated by ADIC /*@ end adic_extract @*/

Applications

We have tested the mixed approach on two applications. The following describes each application briefly and provides the results obtained by using the mixed approach.

Periodic adsorption process: The periodic adsorption process was studied from an optimization point of view in [?,?]. A system of PDAEs in time and space with periodic boundary conditions models the cyclic steady state of a process, where a fluid is preferentially absorbed on the surface of a sorbent bed. This leads to dense Jacobians that dominate the computation time (see [?]). Therefore, previous works have used inexact Jacobians (for example, [?]). Using AD, however, we compute the equality and inequality constraint Jacobians as well as the objective gradient exactly up to machine precision.

The PDAE system is discretized in space by using a finite-volume approach, and the resulting system of ODEs is then integrated in time by using a Runge-Kutta method. This Runge-Kutta iteration in the implementation was determined to be a suitable EDF for differentiation by ADIC, particularly because this routine is repeatedly called at each time step of the simulation and has a C-like implementation. The annotations for declaring this interface routine are shown in Fig. 3(b). Two other lower-level routines for computing the righthand side of the ODE system are also processed by ADIC. The overall problem size depends on the spatial and temporal discretization (N space and N time ). In Table 1 we show the memory required by the trace files created on disk in a purely ADOL-C implementation for various problem sizes, which are absent in the mixed approach. Both approaches preallocate memory of size 2.3 GB in all cases. The absence of trace files on disk in the mixed approach shows that the trace was small enough in all cases to fit into the preallocated memory. Additionally, the runtimes of the mixed approach show improvement over a pure ADOL-C implementation. In Fig. 4(a) the runtimes required in the computation of a equality constraint Jacobian with forward mode are plotted in the left figure for certain problem sizes. In the right side is the time required in the mixed approach is divided into the time required in the wrapper code of the EDF and the ADIC-processed part of the EDF. The same runtimes for the computation of the inequality constraint Jacobian and the objective gradient are shown in Fig. 4(b).

Fluid dynamics -airfoil simulation: Recently, AD was successfully applied to the open source multiphysics suite SU2 [?], which uses a highly modular C++ code structure, to design an efficient adjoint solver [?] for optimization. The implementation is based on the fixed-point formulation of the underlying solver and requires only the recording of one iteration using the converged flow solution. Therefore, at least for steady-state problems, the necessity for checkpointing is eliminated. Still, because of the nature of operator overloading, the memory requirements increase by approximately a factor of 10 compared with the direct flow solver.

SU2 is based on a finite-volume method and offers several well-established combinations of spatial and temporal methods for discretizing the flow equations. Either the steady Euler or the Navier-Stokes equation can be used as the physical model. For this work we have used a second-order central discretization plus an artificial dissipation term (Jameson-Schmidt-Turkel scheme, JST) for the convective terms and a least-squares method for evaluating the gradients needed for the viscous terms. The explicit Euler method is used to advance in pseudotime until convergence. The following two routines were identified as promising use cases for the mixed approach:

1. CCentJST_Flow::ComputeResidual(su2double*val_residual):

per edge convective residual, projects convective flux on the cell-face normal.

2. CEulerSolver::SetPrimitive_Gradient_LS(CGeometry *geometry):

per node gradient of non-conservative variables using least-squares (only Navier-Stokes).

Both routines contain mainly C-like code, which can be processed by ADIC. A potential drawback, however, is that they use mainly class member variables as input. Another difficulty is posed by calls of routines that return variables from other class objects. In such cases we manually copy the data back and forth into simple arrays and define interface routines that take extra inputs.

Figure 5(a) shows the runtime and memory requirements for the Euler solver with a 2D airfoil in transonic flow with 10,216 elements. While the time for tracing is clearly reduced, the evaluation time has significantly increased. This indicates that ADIC-generated derivative code is slower for that case compared with ADOL-C. However there is a decrease of disk usage, solely due to trace files, and an increase in the used RAM. For a 3D airfoil with 582,752 elements the workload for each element is much higher. In that case the total runtime for tracing plus evaluation decreases in the mixed approach as shown in Fig. 5(b). Still, a large fraction of the evaluation time is in the ADIC-generated code. Disk usage reduces by approximately 10% while the used RAM increases insignificantly. For the Navier-Stokes solver with a 2D airfoil with 13937 elements, as shown in Fig. 5(c), the time for tracing reduces by 36%, and the evaluation time increases slightly, resulting in a total reduction by 16%. Furthermore, the total memory usage decreases by 15%.

Conclusions and future work

We have implemented a mixed approach to AD that uses the operator overloading approach to differentiate most of an application and source transformation to differentiate just the computationally intensive portions. The user identifies these portions to be processed by ADIC and uses annotations and a preprocessor to generate code that interfaces ADOL-C's internal data structures with ADIC generated code. The mixed approach has been applied successfully to mediumsized and large-sized applications, resulting in lower memory usage. We plan to apply the mixed approach to more applications. We will also study the benefit of differentiating most of an application using source transformation and only some C++ portions using ADOL-C.

Fig. 1 .

 1 Fig. 1. ADOL-C trace of a simple and externally differentiated function (left). Repeated evaluation of an external function in forward and reverse mode (right)

Fig. 2 .

 2 Fig. 2. (a) Signatures of the forward-and reverse-mode wrapper routines; (b) perroutine registration of EDF; (c) calling an EDF in ADOL-C instrumented code

Fig. 3 .

 3 Fig. 3. (a) Annotations for extracting code for ADIC processing; (b) annotations describing the interface routine to generate wrapper code

Fig. 4 .

 4 Fig. 4. Time required (in seconds) to compute (a) Jacobian of equality constraints in forward mode; (b) Jacobian of inequality constraints and gradient of objective in reverse mode

Fig. 5 .

 5 Fig. 5. Runtime and memory requirements for a 2D Euler case (a), 3D Euler case (b) and 2D Navier-Stokes case (c).

  // primal function signature int myfunc_v2 (int iArrLen, int *iArr, int nin, int nout, int *insz, double **x, int *outsz, double **y, void* ctx); // first order forward implemetation signature int myfunc_forward_v2(int iArrLen, int* iArr, int nin, int nout, int *insz, double **x, int ndir, double ***Xp, int *outsz, double **y, double ***Yp, void* ctx); // first order reverse implementation signature int myfunc_reverse_v2(int iArrLen, int* iArr, int nout, int nin, int *outsz, int dir, double ***Up, int *insz, double ***Zp, double **x, double **y, void* ctx);

		(a)
	ext_diff_fct_v2 * reg_ext_fct_myfunc(){
	ext_diff_fct_v2 *edf	trace_on(tag);
	= reg_ext_fct(myfunc_v2);	... // evaluations
	edf->zos_forward = myfunc_v2;	if (firsttime) edf = reg_ext_fct_myfunc();
	edf->fov_forward = myfunc_forward_v2;	call_ext_fct(edf,...);
	edf->fov_reverse = myfunc_reverse_v2;	... // further evaluations
	// similar for scalar modes	trace_off();
	...	...
	return edf;	jacobian(tag,...); // when required
	}
		(c)
	(b)

Table 1 .

 1 Sizes of trace files created on disk in a purely ADOL-C implementation of periodic adsorbtion process that are not present in a mixed approach

	Nspace Ntime	2000	3000	5000
	20	576 MB	863 MB	2511 MB
	30	856 MB	2240 MB	3734 MB
	50	2472 MB	3707 MB	7146 MB

Acknowledgments. This work was funded in part by a grant from DAAD Project Based Personnel Exchange Programme and by a grant from the U.S. Department of Energy, Office of Science, under contract DE-AC02-06CH11357.

The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory ("Argonne"). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up, nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government.


Martin Benning

Florian Knoll
email: florian.knoll@nyumc.org

Carola-Bibiane Schönlieb

Tuomo Valkonen
email: tuomo.valkonen@liverpool.ac.uk

Preconditioned ADMM with nonlinear operator constraint

Keywords: ADMM, primal-dual, nonlinear inverse problems, parallel MRI, proximal point method, operator splitting, iterative Bregman method

We are presenting a modification of the well-known Alternating Direction Method of Multipliers (ADMM) algorithm with additional preconditioning that aims at solving convex optimisation problems with nonlinear operator constraints. Connections to the recently developed Nonlinear Primal-Dual Hybrid Gradient Method (NL-PDHGM) are presented, and the algorithm is demonstrated to handle the nonlinear inverse problem of parallel Magnetic Resonance Imaging (MRI).

Introduction

Non-smooth regularisation methods are popular tools in the imaging sciences. They allow to promote sparsity of inverse problem solutions with respect to specific representations; they can implicitly restrict the null-space of the forward operator while guaranteeing noise suppression at the same time. The most prominent representatives of this class are total variation regularisation [START_REF] Leonid I Rudin | Nonlinear total variation based noise removal algorithms[END_REF] and 1 -norm regularisation as in the broader context of compressed sensing [START_REF] David | Compressed sensing. Information Theory[END_REF][START_REF] Emmanuel | Compressive sampling[END_REF].

In order to solve convex, non-smooth regularisation methods with linear operator constraints computationally, first-order operator splitting methods have gained increasing interest over the last decade, see [START_REF] Gabay | Applications of the method of multipliers to variational inequalities[END_REF][START_REF] Goldstein | The split bregman method for l1-regularized problems[END_REF][START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] to name just a few. Despite some recent extensions to certain types of non-convex problems [START_REF] Ochs | ipiano: inertial proximal algorithm for nonconvex optimization[END_REF][START_REF] Möllenhoff | The primal-dual hybrid gradient method for semiconvex splittings[END_REF][START_REF] Bonettini | Variable metric inexact line-search based methods for nonsmooth optimization[END_REF][START_REF] Möller | Variational depth from focus reconstruction[END_REF] there has to our knowledge only been made little progress for nonlinear operators constraints [START_REF] Bachmayr | Iterative total variation schemes for nonlinear inverse problems[END_REF][START_REF] Valkonen | A primal-dual hybrid gradient method for nonlinear operators with applications to mri[END_REF].

In this paper we are particularly interested in minimising non-smooth, convex functionals with nonlinear operator constraints. This model covers many interesting applications; one particular application that we are going to address is the joint reconstruction of the spin-proton-density and coil sensitivity maps in parallel MRI [START_REF] Uecker | Image reconstruction by regularized nonlinear inversionjoint estimation of coil sensitivities and image content[END_REF][START_REF] Knoll | Parallel imaging with nonlinear reconstruction using variational penalties[END_REF].

The paper is structured as follows: we will introduce the generic problem formulation, then address its numerical minimisation via a generalised ADMM method with linearised operator constraints. Subsequently we will show connections to the recently proposed NL-PDHGM method (indicating a local convergence result of the proposed algorithm) and conclude with the joint spin-protondensity and coil sensitivity map estimation as a numerical example.

Problem formulation

We consider the following generic constrained minimisation problem:

(û, v) = arg min u,v {H(u) + J(v) subject to F (u, v) = c} . ( 1
)
Here H and J denote proper, convex and lower semi-continuous functionals, F is a nonlinear operator and c a given function. Note that for nonlinear operators of the form F (u, v) = G(u) -v and c = 0 problem (1) can be written as

û = arg min u {H(u) + J(G(u))} . ( 2
)
In the following we want to propose a strategy for solving (1) that is based on simultaneous linearisation of the nonlinear operator constraint and the solution of an inexact ADMM problem.

Alternating direction method of multipliers

We solve (1) by alternating optimisation of the augmented Lagrange function

L δ (u, v; µ) = H(u) + J(v) + µ, F (u, v) -c + δ 2 F (u, v) -c 2 2 . ( 3
)
Alternating minimisation of (3) in u, v and subsequent maximisation of µ via a step of gradient ascent yields this nonlinear version of ADMM [START_REF] Gabay | Applications of the method of multipliers to variational inequalities[END_REF]:

u k+1 ∈ arg min u δ 2 F (u, v k ) -c 2 2 + µ k , F (u, v k ) + H(u) , (4)
v k+1 ∈ arg min v δ 2 F (u k+1 , v) -c 2 2 + µ k , F (u k+1 , v) + J(v) , (5)
µ k+1 = µ k + δ F (u k+1 , v k+1 ) -c . (6)
Not having to deal with nonlinear subproblems, we replace F (u, v k ) and F (u k+1 , v) by their Taylor linearisations around u k and v k , which yields

F (u, v k ) ≈ F (u k , v k )+ ∂ u F (u k , v k ) u -u k and F (u k+1 , v) ≈ F (u k+1 , v k )+∂ v F (u k+1 , v k ) v -v k , re- spectively.
The updates (4) and ( 5) modify to

u k+1 ∈ arg min u δ 2 A k u -c k 1 2 2 + µ k , A k u + H(u) , (7)
v k+1 ∈ arg min v δ 2 B k v -c k 2 2 2 + µ k , B k v + J(v) , (8)
with

A k := ∂ u F (u k , v k ), B k := ∂ v F (u k+1 , v k ), c k 1 := c + A k u k -F (u k , v k ) and c k 2 := c+B k v k -F (u k+1 , v k ).
Note that the updates [START_REF] Bonettini | Variable metric inexact line-search based methods for nonsmooth optimization[END_REF] and ( 8) are still implicit, regardless of H and J. In the following, we want to modify the updates such that they become simple proximity operations.

Preconditioned ADMM

Based on [START_REF] Zhang | A unified primal-dual algorithm framework based on bregman iteration[END_REF], we modify [START_REF] Bonettini | Variable metric inexact line-search based methods for nonsmooth optimization[END_REF] and ( 8) by adding the surrogate terms u k+1 -

u k 2 Q k 1 /2 and v k+1 -v k 2 Q k 2 /2, with w Q :=
Qw, w (note that if Q is chosen to be positive definite, • Q becomes a norm). We then obtain

u k+1 ∈ arg min u δ 2 A k u -c k 1 2 2 + µ k , A k u + H(u) + 1 2 u -u k 2 Q k 1 , v k+1 ∈ arg min v δ 2 B k v -c k 2 2 2 + µ k , B k v + J(v) + 1 2 v -v k 2 Q k 2 . If we choose Q k 1 := τ k 1 I -δA k * A k with τ k 1 δ < 1/ A k 2 and Q k 2 := τ k 2 I -δB k * B k with τ k 2 δ < 1/ B k 2 and if we define µ k := 2µ k -µ k-1 we obtain u k+1 = I + τ k 1 ∂H -1 u k -τ k 1 A k * µ k , (9)
v k+1 = I + τ k 2 ∂J -1 v k -τ k 2 B k * µ k + δ F (u k+1 , v k ) -c , (10)
with (I + α∂E) -1 (w) denoting the proximity or resolvent operator

(I + α∂E) -1 (w) := arg min u 1 2 u -w 2 2 + αE(u) .
The entire proposed algorithm with updates ( 9), ( 10) and ( 6) reads as

Algorithm 1 Preconditioned ADMM with nonlinear operator constraint

Parameters: H, J, F, c Initialization: u 0 , v 0 , µ 0 , δ µ 0 = µ 0 while convergence criterion is not met do

A k = ∂uF (u k , v k ) Set τ k 1 such that τ k 1 δ < 1/ A k 2 u k+1 = I + τ k 1 ∂H -1 u k -τ k 1 A k * µ k B k = ∂vF (u k+1 , v k ) Set τ k 2 such that τ k 2 δ < 1/ B k 2 v k+1 = I + τ k 2 ∂J -1 v k -τ k 2 B k * µ k + δ F (u k+1 , v k ) -c µ k+1 = µ k + δ F (u k+1 , v k+1 ) -c µ k+1 = 2µ k+1 -µ k end while return u k , v k , µ k , µ k 5 Connection to NL-PDHGM
In the following we want to show how the algorithm simplifies in case the nonlinear operator constraint is only nonlinear in one variable, which is sufficient for problems of the form [START_REF] Bachmayr | Iterative total variation schemes for nonlinear inverse problems[END_REF]. Without loss of generality we consider constraints of the form F (u, v) = G(u) -v, where G represents a nonlinear operator in u. Then we have A k = J G(u k ) (with J G(u k ) denoting the Jacobi matrix of G at u k ), B k = -I and if we further choose τ k 2 = 1/δ for all k, update (10) reads

v k+1 = I + 1 δ ∂J -1 G(u k+1 ) + 1 δ µ k . Applying Moreau's identity [18] b = I + 1 δ ∂J -1 (b) + 1 δ (I + δ∂J * ) -1 (δb) yields µ k+1 = (I + δ∂J * ) -1 µ k + δG(u k+1 ) .
If we further change the order of the updates, starting with the update for µ, the whole algorithm reads

µ k+1 = (I + δ∂J * ) -1 µ k + δG(u k ) , µ k+1 = 2µ k+1 -µ k , u k+1 = I + τ k 1 ∂H -1 u k -τ k 1 J G(u k ) * µ k+1 .
Note that this algorithm is almost the same as NL-PDHGM proposed in [START_REF] Valkonen | A primal-dual hybrid gradient method for nonlinear operators with applications to mri[END_REF] for θ = 1, except that the extrapolation step is carried out on the dual variable µ instead of the primal variable u. In the following we want to briefly sketch how to prove convergence for this algorithm in analogy to [START_REF] Valkonen | A primal-dual hybrid gradient method for nonlinear operators with applications to mri[END_REF]. We define

N (µ k+1 , u k+1 ) := ∂J * (µ k+1 ) -∇G(u k )u k+1 -c k ∂H(u k+1 ) + J G(u k ) * µ k+1 , L k := 1 δ I J G(u k ) J G(u k ) * 1 τ k 1 I , with c k := G(u k ) -J G(u k )u k . Now the algorithm is: find (µ k+1 , u k+1 ) such that N (µ k+1 , u k+1 ) + L k (µ k+1 -µ k , u k+1 -u k ) 0.
If we exchange the order of µ and u here, i.e., reorder the rows of N , and the rows and columns of L k , we obtain almost the "linearised" NL-PDHGM of [START_REF] Valkonen | A primal-dual hybrid gradient method for nonlinear operators with applications to mri[END_REF]. The difference is that the sign of J G in L k is inverted. The only points in [START_REF] Valkonen | A primal-dual hybrid gradient method for nonlinear operators with applications to mri[END_REF] where the exact structure of L k (M x k therein) is used, are Lemma 3.1, Lemma 3.6 and Lemma 3.10. The first two go through exactly as before with the negated structure. Reproducing Lemma 3.10 demands bounding actual step lengths u k -u k+1 and µ k -µ k+1 from below, near a solution for arbitrary > 0. A proof would go beyond the page limit of this proceeding. Let us just point out that this can be done, implying that the convergence results of [START_REF] Valkonen | A primal-dual hybrid gradient method for nonlinear operators with applications to mri[END_REF] apply for this algorithm as well. This means that under somewhat technical regularity conditions, which for TV type problems amount to Huber regularisation, local convergence in a neighbourhood of the true solution can be guaranteed.

Joint estimation of the spin-proton density and coil sensitivities in parallel MRI

We want to demonstrate the numerical capabilities of Algorithm 1 by applying it to the nonlinear problem of joint estimation of the spin-proton density and the coil sensitivities in parallel MRI. The discrete problem of joint reconstruction from sub-sampled k-space data on a rectangular grid reads

     û ĉ1 . . . ĉ2      ∈ arg min v=(u,c1,...,cn)    1 2 n j=1 SF(G(v)) j -f j 2 2 + α 0 R 0 (u) + n j=1 α j R j (c j )    ,
where F is the 2D discrete Fourier transform, f j are the k-space measurements for each of the n coils, S is the sub-sampling operator and R j denote appropriate regularisation functionals. The nonlinear operator G maps the unknown spinproton density u and the different coil sensitivities c j as follows [START_REF] Uecker | Image reconstruction by regularized nonlinear inversionjoint estimation of coil sensitivities and image content[END_REF]:

G(u, c 1 , . . . , c n ) = (uc 1 , uc 2 , . . . , uc n ) T . (11)
In order to compensate for sub-sampling artefacts in sub-sampled MRI it is common practice to use total variation as a regulariser [START_REF] Block | Undersampled radial mri with multiple coils. iterative image reconstruction using a total variation constraint[END_REF][START_REF] Ramani | Parallel mr image reconstruction using augmented lagrangian methods[END_REF]. Coil sensitivities are assumed to be smooth, cf. Figure 1, motivating a reconstruction model similar to the one proposed in [START_REF] Knoll | Parallel imaging with nonlinear reconstruction using variational penalties[END_REF]. We therefore choose the discrete isotropic total variation, R 0 (u) = ∇u 2,1 , and the smooth 2-norm of the discretised gradient, i.e. R j (c j ) := ∇c j 2,2 , for all j > 0, following the notation in [START_REF] Benning | Phase reconstruction from velocity-encoded mri measurementsa survey of sparsity-promoting variational approaches[END_REF]. We further introduce regularisation parameters λ j in front of the data fidelities and rescale all regularisation parameters such that α 0 + 1 n n j=1 λ j + n j=1 α j = 1. In order to realise this model via Algorithm 1 we consider the following operator splitting strategy. We define F (u 0 , . . . , u n , v 0 , . . . , v 2n ) as

F (u 0 , . . . , u n , v 1 , . . . , v n ) :=         G(u 0 , . . . , u n ) ∇u 0 0 • • • 0 0 ∇u 1 . . . . . . . . . . . . . . . 0 0 • • • 0 ∇u n         -         v 0 . . . v n . . . v 2n        
, set H(u 0 , . . . , u n ) ≡ 0, and

J(v 0 , . . . , v 2n ) = 2n j=0 J j (v j ) with J j (v j ) := λj 2 SFv j - f j 2 2 for j ∈ {0, . . . , n -1}, J n (v n ) = α 0 v n 2,
1 and J j (v j ) = α j-n v j 2,2 for j ∈ {n + 1, . . . , 2n}. Note that with these choices of functions, all the resolvent operations can be carried out easily. In particular, we obtain

(I + τ k 1 ∂H) -1 (w) = w , (I + τ k 2 ∂J j ) -1 (w) = F -1 Fw j + τ k 2 λ j S T f j 1 + τ k 2 λ j diag(S T S)
for j ∈ {0, . . . , n -1} ,

(I + τ k 2 ∂J n ) -1 (w) = w n w n 2 max w n 2 -α 0 τ k 2 , 0 , (I + τ k 2 ∂J j ) -1 (w) = w j w j 2,2 max w j 2,2 -α j-n τ k 2 , 0 for j ∈ {n + 1, . . . , 2n} .
Moreover, as B k = -I (and thus, B k = 1) for all k, we can simply eliminate τ k 2 by replacing it with 1/δ, similar to Section 5. Fig. 1: Figure 1a shows the brain phantom as described in Section 6.1. Figure 1c 1j show visualisations of the measured coil sensitivities of a water bottle. Figure 1b shows the simulated, spiral-shaped sub-sampling scheme used to sub-sample the k-space data. Fig. 2: Reconstructions for noise with low noise level σ = 0.05. Despite the subsampling, features of the brain phantom are very well preserved. In addition, the coil sensitivities seem to correspond well to the original ones, despite a slight loss of contrast. Note that coil sensitivities remain the initial value where the signal is zero.

Experimental setup

We now want to discuss the experimental setup. We want to reconstruct the synthetic brain phantom in Figure 1a from sub-sampled k-space measurements. The numerical phantom is based on the design in [START_REF] Aubert-Broche | A new improved version of the realistic digital brain phantom[END_REF] with a matrix size of 190 × 190. It consists of several different tissue types like cerebrospinal fluid (CSF), gray matter (GM), white matter (WM) and cortical bone. Each pixel is assigned a set of MR tissue properties: Relaxation times T 1 (x, y) and T 2 (x, y) and spin density ρ(x, y). These parameters were also selected according to [START_REF] Aubert-Broche | A new improved version of the realistic digital brain phantom[END_REF]. The MR signal s(x, y) in each pixel was then calculated by using the signal equation of a fluid attenuation inversion recovery (FLAIR) sequence [START_REF] Bernstein | Handbook of MRI pulse sequences[END_REF]:

s(x, y) = ρ(x, y)(1 -2 e -TI/T1(x,y) )(1 -e -TR/T1(x,y) ) e -TE/T2(x,y) .

The sequence parameters were selected: TR = 10000 ms, TE = 90 ms. TI was set to 1781 ms to achieve signal nulling of CSF (T csf 1 log(2) with T csf 1 = 2569ms). In order to generate artificial k-space measurements for each coil, we proceed as follows. First, we produce 8 images of the brain phantom multiplied by the measured coil sensitivity maps shown in Figure 1c -1j. The coil sensitivity maps were generated from the measurements of a water bottle with an 8-channel head coil array. Then we produce artificial k-space data by applying the 2D discrete Fourier-transform to each of those individual images. Subsequently, we sub-sample only approx. 25% of each of the k-space datasets via the spiral shown in Figure 1b. Finally, we add Gaußian noise with standard deviation σ to the sub-sampled data.

Computations

For the actual computations we use two noisy versions f j of the simulated kspace data; one with small noise (σ = 0.05) and one with a high amount of noise (σ = 0.95). As stopping criterion we simply choose a fixed number of iterations; for both the low noise level as well as the high noise level dataset we have fixed the number of iterations to 1500. The initial values used for the algorithm are u 0 j = 1 with 1 ∈ R l×1 being the constant one-vector, for all j ∈ {0, . . . , n}. All other initial variables (v 0 , µ 0 , µ 0 ) are set to zero.

Low noise level

We have computed reconstructions from the noisy data with noise level σ = 0.05 via Algorithm 1, with regularisation parameters set to λ j = 0.0621, α 0 = 0.0062 and α j = 0.9317 for j ∈ {1, . . . , n}. We have further created a naïve reconstruction by averaging the individual inverse Fourier-transformed images obtained from zero-filling the k-space data. The modulus images of the results are visualised in Figure 2. The PSNR values for the averaged zero-filled reconstruction is 10.2185, whereas the PSNR of the reconstruction with the proposed method is 24.5572.

High noise level We proceeded as in the previous section, but for noisy data with noise level σ = 0.95. The regularisation parameters were set to λ j = 0.0149, α 0 = 0.0135 and α j = 0.9716 for j ∈ {1, . . . , n}. The modulus images of the results are visualised in Figure 3. The PSNR values for the averaged zero-filled reconstruction is 9.9621, whereas the PSNR of the reconstruction with the proposed method is 16.672.

Conclusions & outlook

We have presented a novel algorithm that allows to compute minimisers of a sum of convex functionals with nonlinear operator constraint. We have shown the connection to the recently proposed NL-PDHGM algorithm which implies local convergence results in analogy to those derived in [START_REF] Valkonen | A primal-dual hybrid gradient method for nonlinear operators with applications to mri[END_REF]. Subsequently we have demonstrated the computational capabilities of the algorithm by applying it to a nonlinear joint reconstruction problem in parallel MRI.

For future work, the convergence of the algorithm in the general setting has to be verified, and possible extensions to guarantee global convergence have to be studied. Generalisation of stopping criteria such as a linearised primal-dual gap will also be of interest as well. With respect to the presented parallel MRI application, exact conditions for the convergence (like the exact norm of the bounds) have to be verified. The impact of the algorithm-as well as the regularisationparameters on the reconstruction has to be analysed, and a rigorous study with artificial and real data would also be desirable. Moreover, future research will focus on alternative regularisation functions, e.g. based on spherical harmonics motivated by [START_REF] Sbrizzi | Robust reconstruction of b1+ maps by projection into a spherical functions space[END_REF]. Last but not least, other applications that can be modelled via (1) should be considered in future research.

Fig. 3 :

 3 Fig.3: Reconstructions for noise with high noise level σ = 0.95. Due to the large amount of noise, higher regularisation parameters are necessary. As a consequence, fine structures are smoothed out and in contrast to the case of little noise, compensation of sub-sampling artefacts is less successful.

Acknowledgments MB, CS and TV acknowledge EPSRC grant EP/M00483X/1. FK ackowledges National Institutes of Health grant NIH P41 EB017183.

EPSRC data statement: the corresponding code and data are available for download at https://www.repository.cam.ac.uk/handle/1810/256221


Andrei Dmitruk
email: dmitruk@member.ams.org

Nikolai Osmolovskii
email: osmolovski@uph.edu.pl

Necessary optimality conditions in a problem with integral equations on a nonfixed time interval subject to mixed and state constraints

Keywords: Volterra-type equation, extended weak minimum, local maximum principle, state-control constraints, adjoint equation, transversality conditions, change of time variable, linear-positive independence

We consider an optimal control problem with Volterra-type integral equations on a nonfixed time interval subject to endpoint constraints, mixed state-control constraints of equality and inequality type, and pure state inequality constraints. The main assumption is the linearpositive independence of the gradients of active mixed constraints with respect to the control. We formulate first order necessary optimality conditions for an extended weak minimum, the notion of which is a natural generalization of the notion of weak minimum with account of variations of the time. The presented conditions generalize the local maximum principle in optimal control problems with ordinary differential equations.

Introduction

The results presented in this paper generalize the results obtained in our previous two papers [START_REF] Dmitruk | Necessary conditions for a weak minimum in optimal control problems with integral equations subject to state and mixed constraints[END_REF] and [START_REF] Dmitruk | Necessary conditions for a weak minimum in optimal control problems with integral equations on a variable time interval[END_REF]. Paper [START_REF] Dmitruk | Necessary conditions for a weak minimum in optimal control problems with integral equations subject to state and mixed constraints[END_REF] was devoted to the first order necessary conditions for a weak minimum in a general optimal control problem with Volterratype integral equations, considered on a fixed time interval, subject to endpoint constraints of equality and inequality type, mixed state-control constraints of inequality and equality type, and pure state constraints of inequality type. Paper [START_REF] Dmitruk | Necessary conditions for a weak minimum in optimal control problems with integral equations on a variable time interval[END_REF] studied first order necessary conditions for an extended weak minimum in an optimal control problem with Volterra-type integral equations considered on a non-fixed time interval, subject to endpoint constraints of equality and inequality type, but without mixed state-control constraints and pure state constraints.

Here we consider a problem generalizing both problems of [START_REF] Dmitruk | Necessary conditions for a weak minimum in optimal control problems with integral equations subject to state and mixed constraints[END_REF] and [START_REF] Dmitruk | Necessary conditions for a weak minimum in optimal control problems with integral equations on a variable time interval[END_REF]. We formulate first order necessary conditions for an extended weak minimum in this general problem. Following the tradition, we call them stationarity conditions, or conditions of the local maximum principle . They are presented in Theorem 1.

As far as we know, such conditions for problems with integral equations on a variable time interval were not obtained up to now. Their novelty, as compared with those for problems on a fixed time interval is that the costate equation and transversality condition with respect to t involve nonstandard terms that are absent in problems with ODEs. More remarks concerning the existing literature on the problems with integral equations can be found in papers [1-4, 6, 7].

As was already mentioned in [START_REF] Dmitruk | Necessary conditions for a weak minimum in optimal control problems with integral equations subject to state and mixed constraints[END_REF], the stationarity conditions in optimal control problems constitute an important stage in obtaining any further necessary optimality condition, including maximum principle or higher order conditions, and thus, they deserve a separate thorough study for each specific class of problems.

The paper is organized as follows. In Section 2 we formulate a general optimal control problem with integral equations on a variable time interval which we call Problem A. We also define in this section the notion of the extended weak minimum. Section 3 is devoted to formulation of the main result of the paper -the local maximum principle in Problem A, which is the first order necessary condition for an extended weak minimum (Theorem 1). A short discussions of its proof is given in Section 4.

General optimal control problem with integral equations on a variable time interval (Problem A)

Consider the following control system of Volterra-type integral equations on a variable time interval [t 0 , t 1 ] :

x(t) = x(t 0 ) + t t0 f (t, s, x(s), u(s)) ds, (1)
where x(•) is a continuous n-dimensional and u(•) is a measurable essentially bounded r-dimensional vector-functions on [t 0 , t 1 ]. As usual, we call x(•) the state variable and u(•) the control variable (or simply the control ). A pair w(t) = (x(t), u(t)) defined on its own interval [t 0 , t 1 ] and satisfying (1) for a.e. t ∈ [t 0 , t 1 ] is called a process. We assume that the function f is defined and twice continuously differentiable on an open set R ⊂ R 2+n+r .

The problem is to minimize the endpoint functional

J = ϕ 0 (t 0 , x(t 0 ), t 1 , x(t 1 )) → min (2)
on the set of all processes (solutions of system (1)) satisfying the endpoint constraints

η j (t 0 , x(t 0 ), t 1 , x(t 1 )) = 0, j = 1, . . . , k, (3)
ϕ i (t 0 , x(t 0 ), t 1 , x(t 1 )) 0, i = 1, . . . , ν, (4)
the mixed state-control constraints

F i (t, x(t), u(t)) 0 for a.e. t ∈ [t 0 , t 1 ], i = 1, . . . , d(F ), (5)
G j (t, x(t), u(t)) = 0 for a.e. t ∈ [t 0 , t 1 ], j = 1, . . . , d(G), (6)
and the state constraints

Φ k (t, x(t)) 0 for all t ∈ [t 0 , t 1 ], k = 1, . . . , d(Φ). (7)
The functions ϕ 0 , ϕ i , η j are assumed to be defined and continuously differentiable on an open set P ⊂ R 2n+2 , and the functions F i , G j , and Φ k are assumed to be defined and continuously differentiable on an open set Q ⊂ R 1+n+r (the smoothness assumptions). The notation d(F ), d(G), and d(Φ) stand for the numbers of these functions. Moreover, we assume that the mixed constraints ( 5) and ( 6) are regular in the following sense: at any point (t, x, u) ∈ Q satisfying relations F i 0 ∀ i and G j = 0 ∀ j, the system of vectors

F iu (t, x, u), i ∈ I(t, x, u), G ju (t, x, u), j = 1, . . . , d(G),
is positively-linearly independent, where I(t, x, u) = {i | F i (t, x, u) = 0 } is the set of active indices of mixed inequality constraints at the given point. Here and in the sequel we denote by F iu the partial derivative (gradient) of the function F i with respect to the variable u. Similar notation is used for other functions an variables. Recall that a system consisting of two tuples of vectors p 1 , . . . , p m and q 1 , . . . q k in the space R r is said to be positively-linearly independent if there does not exist a nontrivial tuple of multipliers α 1 , . . . , α m , β 1 , . . . β k with all α i 0 such that

i α i p i + j β j q j = 0.
The problem (1)-( 7) will be called Problem A, and the relations (2)-( 4) its endpoint block.

Note that the function f explicitly depends on two time variables, t and s, the roles of which are essentially different. Conventionally, the variable s will be called inner, while t will be called outer time variable, and one should carefully distinguish between them in further considerations. Among the four arguments of the function f and its derivatives, the first argument will always be the outer and the second one be the inner time variables, no matter by which letters they will be denoted.

As in [START_REF] Dmitruk | Necessary conditions for a weak minimum in optimal control problems with integral equations subject to state and mixed constraints[END_REF] and [START_REF] Dmitruk | Necessary conditions for a weak minimum in optimal control problems with integral equations on a variable time interval[END_REF], we mention an important particular case of system (1): if f does not depend on the outer time variable t, i.e., f = f (s, x(s), u(s)), then the integral equation ( 1) is equivalent to the differential equation ẋ(t) = f (t, x(t), u(t)), hence Problem A becomes an optimal control problem of ordinary differential equations on a nonfixed time interval.

Obviously, each process under consideration must "lie" in the domain R of the function f (t, s, x, u), i.e.

(t, s, x(s), u(s)) ∈ R for a.e. (t, s) ∈ ∆[t 0 , t 1 ],
where ∆[t 0 , t 1 ] = { (t, s) : t 0 s t t 1 }. We will need even a stronger condition.

Definition. A process w(t) = (x(t), u(t)) defined on an interval t ∈ [t 0 , t 1 ] (with continuous x(t) and measurable and essentially bounded u(t) ) will be called admissible with respect to R if its "extended graph"

Γ (w) = {(t, s, x(s), u(s)) | (t, s) ∈ ∆[t 0 , t 1 ] } lies in the set R with some "margin" , i.e., dist ((t, s, x(s), u(s)), ∂R) const > 0 for a.a. (t, s) ∈ ∆[t 0 , t 1 ].
A process is called admissible in problem A if it is admissible with respect to R and satisfies all the constraints (1) and ( 3)-( 7) of the problem.

Like in any problem on a nonfixed time interval, the notion of weak minimum in Problem A needs a modification.

Definition. We will say that an admissible process w 0 (t) = (x 0 (t), u 0 (t)), t ∈ [ t0 , t1 ], provides the extended weak minimum if there exists an ε > 0 such that for any Lipschitz continuous bijective mapping ρ : [t 0 , t 1 ] → [ t0 , t1 ] satisfying the conditions |ρ(t) -t| < ε and | ρ(t) -1| < ε, and for any admissible process w(t) = (x(t), u(t)), t ∈ [t 0 , t 1 ], satisfying the conditions

|x(t) -x 0 (ρ(t))| ε ∀ t, and |u(t) -u 0 (ρ(t))| ε (∀ )t, (8)
we have J(w) J(w 0 ). (Notation (∀ ), as usual, means "for almost all".)

The conditions on ρ imply ρ(t 0 ) = t0 and ρ(t 1 ) = t1 with | t0 -t 0 | < ε and | t1 -t 1 | < ε. If the interval [t 0 , t 1 ] is fixed and we take ρ(t) = t, then relations (8) describe the usual uniform closeness between the processes w 0 and w both in the state and control variables. However, for an arbitrary ρ(t), relations (8) extend the set of "competing" processes, and thus, even for a fixed time interval, the extended weak minimum is stronger than the usual weak minimum.

Local Maximum Principle in Problem A

Let a process

w 0 (t) = (x 0 (t), u 0 (t)), t ∈ [ t0 , t1 ] (9)
provide an extended weak minimum in Problem A. We assume that the endpoints of the reference state x 0 (t) do not lie on the boundary of state constraints; moreover, that

Φ k ( t0 , x 0 ( t0 )) < 0, Φ k ( t1 , x 0 ( t1 )) < 0, k = 1, . . . , d(Φ). (10)
For process (9), let us introduce a tuple of Lagrange multipliers

λ = (α 0 , α, β, ψ x (t), ψ t (t), h i (t), m j (t), µ k (t)). ( 11
)
Here α 0 ∈ R corresponds to the cost ϕ 0 , the components of vectors α = (α 1 , . . . , α ν ) ∈ R ν and β = (β 1 , . . . , β k ) ∈ R k correspond to endpoint constraints ϕ i 0, and η j = 0, respectively, adjoint variables ψ x (t) : [ t0 , t1 ] → R n and ψ t (t) : [ t0 , t1 ] → R correspond to the control system, µ k (t) : [ t0 , t1 ] → R, k = 1, . . . , d(Φ), refer to the state constraints Φ k (t, x) 0; moreover, ψ x , ψ t , and µ k are functions of bounded variation, continuous at t0 and t1 ; multipliers h i (t) : [ t0 , t1 ] → R, i = 1, . . . , d(F ), and m j (t) : [ t0 , t1 ] → R, j = 1, . . . , d(G) corresponding to the mixed constraints F i (t, x, u) 0, G j (t, x, u) = 0, are measurable bounded functions.

Note that here ψ x and ψ t are not the partial derivatives with respect to x and t, but simply the adjoint variables, which refer to x and t, respectively. This notation was proposed by Dubovitskii and Milyutin and turned out to be highly convenient, especially in problems with many state variables. We hope it will not cause confusion.

We denote by dψ x , dψ t , dµ k the Lebesgue-Stieltjes measures which correspond to the functions of bounded variation ψ x , ψ t , µ k , respectively. These measures have no atoms at the points t0 and t1 , and moreover, dµ k 0, k = 1, . . . , d(Φ), since it corresponds to the inequality constraint. Hence each µ k is a monotone nondecreasing function. By

dψ x dt = ψx (t), dψ t dt = ψt (t), dµ k dt = μk (t)
we denote the generalized derivatives of these functions with respect to t. Consequently, the following relations hold:

ψx (t) dt = dψ x (t), ψt (t) dt = dψ t (t), μk (t) dt = dµ k (t).
In what follows, all pointwise relations involving continuous functions hold for all t, while those involving measurable functions hold for almost all t.

In order to present optimality conditions in Problem A, introduce, for a tuple λ of (11), the modified Pontryagin function

H(t, s, x, u) = ψ x (t-) f (t, s, x, u) + t1 t ψ x (τ ) f t (τ, s, x, u) dτ. (12)
Here ψ x (t-) means the left hand value of the function ψ x at a point t, and f t means the partial derivative of the function f (t, s, x, u) with respect to the first, outer variable t.

Also, for w 0 and λ, let us introduce the augmented modified Pontryagin function H(t, s, x, u) = H(t, s, x, u)

- i h i (t) F i (s, x, u) - j m j (t) G j (s, x, u) - k μk (t) Φ k (s, x), (13)
the endpoint Lagrange function

l(t 0 , x 0 , t 1 , x 1 ) = ν i=0 α i ϕ i + k j=1
β j η j (t 0 , x 0 , t 1 , x 1 ), and a special auxiliary function

R(t) = t t0 f t (t, s, x 0 (s), u 0 (s)) ds.
The main result of the paper is the following Theorem 1 (local maximum principle). If a process w 0 (t) = (x 0 (t), u 0 (t)), t ∈ [ t0 , t1 ] provides the extended weak minimum in Problem A and satisfies assumption (10), then there exists a tuple of multipliers (11) satisfying the specified above properties and such that the following conditions hold true: a) nonnegativity conditions

α 0 0, α 0, h i (t) 0, i = 1, . . . , d(F ), dµ k 0, k = 1, . . . , d(Φ), b) nontrivality condition α 0 + |α| + |β| + k µ k ( t1 ) -µ k ( t0 ) + i t1 t0 h i (t) dt > 0,
c) endpoint complementary slackness conditions α i ϕ i ( t0 , x 0 ( t0 ), t1 , x 0 ( t1 )) = 0, i = 1, . . . , ν, d) pointwise complementary slackness conditions

dµ k (t) Φ k (t, x 0 (t)) ≡ 0, k = 1, . . . , d(Φ), h i (t) F i (t, x 0 (t), u 0 (t)) = 0 a.e. on [ t0 , t1 ], e) adjoint equation in x -dψ x (t) = H x (t, t, x 0 (t), u 0 (t))dt = ψ x (t-)f x (t, t, x 0 (t), u 0 (t)) + t1 t ψ x (τ ) f tx (τ, t, x 0 (t), u 0 (t)) dτ dt - i h i (t) F ix (t, x 0 (t), u 0 (t)) + j m j (t) G jx (t, x 0 (t), u 0 (t)) dt - k dµ k (t) Φ kx (t, x 0 (t)), f ) adjoint equation in t -dψ t (t) = H s (t, t, x 0 (t), u 0 (t)) dt - i h i (t) F it (t, x 0 (t), u 0 (t)) + j m j (t) G jt (t, x 0 (t), u 0 (t)) dt - k dµ k (t) Φ kt (t, x 0 (t)) -dψ x (t)R(t), ( 14
) with H s (t, t, x 0 (t), u 0 (t)) = ψ x (t-)f s (t, t, x 0 (t), u 0 (t)) + t1 t ψ x (τ ) f ts (τ, t, x 0 (t), u 0 (t)) dτ,
where f s is the partial derivative of the function f (t, s, x, u) with respect to the second, inner variable s, and f ts is its second partial derivative, g) transversality conditions in x,

ψ x ( t0 ) = l x0 , -ψ x ( t1 ) = l x1 , h) transversality conditions in t, ψ t ( t0 ) = l t0 , -ψ t ( t1 ) = l t1 -ψ x ( t1 )R( t1 ), (15)
i) stationarity condition with respect to the control

H u (t, t, x 0 (t), u 0 (t)) = 0 a.e. on [ t0 , t1 ],
i.e.,

ψ x (t-)f u (t, t, x 0 (t), u 0 (t)) + t1 t ψ x (τ ) f tu (τ, t, x 0 (t), u 0 (t)) dτ - i h i (t) F iu (t, x 0 (t), u 0 (t)) + j m j (t) G ju (t, x 0 (t), u 0 (t)) = 0, k) the "energy evolution law" H(t, t, x 0 (t), u 0 (t)) + ψ t (t) = 0 a.e. on [ t0 , t1 ].
The last condition is called in such a way, since together with (14) it gives the equation for evolution of the function H(t, t, x 0 (t), u 0 (t)), which is often (especially in mechanical problems) regarded as the total energy of the system:

Ḣ = H s -ψx R .
If the state and mixed constraints are absent and the dynamics does not explicitly depend on time: f = f (x, u), then H = H, R = 0, and we get "the energy conservation law": H = const along the optimal process. Note that both the adjoint equation in t and the right transversality condition in t involve additional terms, dψ x (t)R(t) and ψ x ( t1 )R( t1 ), respectively. Both these terms are generated by the dependence of f (t, s, x, u) on the outer time variable t, which was absent in problems with ODEs. (Indeed, in those problems f t = 0, whence F = 0, so this term disappears.) In our opinion, this novelty in optimality conditions for problems on a variable time interval needs further study.

Using generalized derivatives of functions of bounded variation, we can represent the adjoint equation in x and t in the easy-to-remember form:

- dψ x (t) dt = H x (t, t, x 0 (t), u 0 (t)) = H x (t, t, x 0 (t), u 0 (t)) - i h i (t) F ix (t, x 0 (t), u 0 (t)) - j m j (t) G jx (t, x 0 (t), u 0 (t)) - k dµ k (t) dt Φ kx (t, x 0 (t)), (16)
and

- dψ t (t) dt = H s (t, t, x 0 (t), u 0 (t)) - i h i (t) F it (t, x 0 (t), u 0 (t)) - j m j (t) G jt (t, x 0 (t), u 0 (t)) - k dµ k (t) dt Φ kt (t, x 0 (t)) - dψ x (t) dt R(t). (17)
4 About the Proof of Theorem 1

Like in our paper [START_REF] Dmitruk | Necessary conditions for a weak minimum in optimal control problems with integral equations on a variable time interval[END_REF], in order to prove Theorem 1, we reduce Problem A to an auxiliary problem on a fixed time interval by using the change of time variable t = t(τ ), where dt/dτ = v(τ ) and v(τ ) > 0. Setting x(τ ) = x(t(τ )) and ũ(τ ) = u(t(τ )), we come to the following system of integral equations:

x(τ ) = x(τ 0 ) + τ τ0 f (t(τ ), t(σ), x(σ), ũ(σ)) v(σ) dσ, (18)
t(τ ) = t(τ 0 ) + τ τ0 v(σ) dσ, ( 19
)
where τ is a new time, t(τ ) an additional state variable, v(τ ) an additional control variable, and σ a new time of integration instead of s.

We see that here the integrand of the first equation involves the value t(τ ) of a state variable, which was not allowed in [START_REF] De La Vega | Necessary conditions for optimal terminal time control problems governed by a Volterra integral equation[END_REF]. Abstracting from the specific form of the second equation and changing the notation t(τ ) to a more general y(τ ) (and changing also τ to a more convenient t ), we come to a system of the following form on a fixed interval [t 0 , t 1 ] : This system does not fall into the framework of equation ( 1), since the integrand of the first equation depends on y(t) (which can be regarded as the outer state variable). Thus, we have to study a new, broader than (1), class of integral control systems.

x(t) = x(t 0 ) + t t0 g(t,
Adding to the obtained system the mixed constraints, the state constraints and the terminal block, we obtain the following

Problem B on a fixed time interval

On the set of solutions w = (x, y, u) to system (20)-( 21) satisfying the constraints

F i (t,
η j (x(t 0 ), y(t 0 ), x(t 1 ), y(t 1 )) = 0, j = 1, . . . , k,

ϕ i (x(t 0 ), y(t 0 ), x(t 1 ), y(t 1 )) 0, i = 1, . . . , ν,

to minimize the endpoint functional J = ϕ 0 (x(t 0 ), y(t 0 ), x(t 1 ), y(t 1 )) → min .

Like before, the functions η j , ϕ i , and ϕ 0 are assumed to be continuously differentiable on an open set P ⊂ R 2n+2m , the functions F i , G j , and Φ k continuously differentiable on an open set Q ⊂ R 1+m+n+r . We also assume that the mixed constraints ( 22)-( 23) are regular in the same sense as in Problem A.

To derive optimality conditions in Problem B, we consider it as a particular case of an abstract nonsmooth problem in a Banach space, hence we can apply the well known abstract Lagrange multipliers rule for nonsmooth problems (see, e.g. [START_REF] Milyutin | Maximum principle in optimal control[END_REF][START_REF] Dmitruk | Necessary conditions for a weak minimum in optimal control problems with integral equations subject to state and mixed constraints[END_REF]). Let us formulate it. We study the local minimality of an admissible point x 0 ∈ D. Assume that the cost f 0 and the mappings b i are Frechet differentiable at x 0 , the operator g is strictly differentiable at x 0 , and the image of g (x 0 ) is closed. Let K 0 i be the polar cone to K i , i = 1, . . . , ν. is stationary at x 0 : L (x 0 ) = 0.

Applying Theorem 2 to Problem B, we perform some analysis of the obtained conditions and represent them in the form of local maximum principle for Problem B. The latter is then applied to the auxiliary problem with system (18)-( 19), and finally, we rewrite the results in terms of the original Problem A.

4. 2

 2 Lagrange multipliers rule for an abstract nonsmooth problem Let X, Y, Z i , i = 1, . . . , ν be Banach spaces, D ⊂ X an open set, K i ⊂ Z i i = 1, . . . , ν closed convex cones with nonempty interiors, f 0 : D → R, b i : D → Z i , i = 1, . . . , ν, and g : D → Y given mappings. Consider the following problem f 0 (x) → min, b i (x) ∈ K i , i = 1, . . . , ν, g(x) = 0. (28)

Theorem 2 .

 2 Let x 0 provide a local minimum in problem (28). Then there exist Lagrange multipliers α 0 0, ζ * i ∈ K 0 i , i = 1, . . . , ν, and y * ∈ Y * , not all equal to zero, satisfying the complementary slackness conditions ζ * i , b i (x 0 ) = 0, i = 1, . . . , ν, and such that the Lagrange functionL(x) = α 0 f 0 (x) + ν i=1 ζ * i , b i (x) + y * , g(x)

  and y(t) are continuous functions of dimensions n and m respectively, u(t) is a measurable and essentially bounded function on [t 0 , t 1 ]. We still denote the time by t. The data functions g and h, as before, are assumed to be twice continuously differentiable on an open set R ⊂ R 2+2m+n+r .

	s, y(t), x(s), y(s), u(s)) ds,	(20)
	t
	y(t) = y(t 0 ) +	h(t, s, y(s), u(s)) ds,	(21)
	t0
	where x(t)