Jacques Palicot
email: jacques.palicot@centralesupelec.fr

Christophe Moy

Benoit Résimont

Rémi Bonnefoi

Application of Hierarchical and Distributed Cognitive Architecture Management for the Smart Grid

Keywords: Smart Grid (SG), Smart home management, Hierarchical and Distributed Cognitive Radio Architecture Management (HDCRAM), multi-agent systems, JADE

Moving from the current power grid to the Smart Grid (SG) requires decentralizing management. This should be done by distributing intelligence over the entire grid, thereby, the intermittent production of renewable energy, customer consumption and electricity storage in electrical vehicles (EVs) could be managed in real time. In this paper, the Hierarchical and Distributed Cognitive Radio Architecture Management (HDCRAM), initially

proposed to manage Cognitive Radio systems, is proposed for the management of the SG. This architecture can both be applied to the whole SG and to any sub-part (distribution network, production network, microgrid). In this paper we focus on the distribution network and the hierarchical position of each element is identified. As an example, HDCRAM is used for smart home management and multi-agent based modeling shows benefits of such an architecture. In the simulated scenario, without any management the peak power consumption is 5500W and the hierarchical and distributed management allows to reduce it to 900 W. This diminution allows to reduce the pressure on the grid and can decrease the risk of failure.

Introduction

The European Technology Platform (ETP) Smart Grids was set up in 2005 to create a joint vision of European networks for 2020 and onwards [START_REF] Commission | European Technology Platform SmartGrids Vision and Strategy for Europes Electricity Networks of the Future[END_REF].

According to the European Commission, the definition of the Smart Grid is a follows [START_REF]Expert Group 1: Functionalities of Smart Grids and Smart Meters: Final Deliverable[END_REF]. Definition 1. A Smart Grid is an electricity network that can cost efficiently integrate the behavior and actions of all users connected to it -generators, consumers and those that do both -in order to ensure economically efficient, sustainable power system with low losses and high levels of quality and security of supply and safety.

A more explicit definition is given by the Energy Independence and Security Act of 2007 [START_REF]110th United States Congress, "Energy Independence and Security Act of 2007[END_REF]. Definition 2. The term "Smart Grid" refers to a modernization of the electricity delivery system so that it monitors, protects, and automatically optimizes the operation of its interconnected elements, from the central and distributed generator through the high-voltage transmission network and the distribution system, to industrial users and building automation systems, to energy storage installations, and to end-use consumers and their thermostats, electric vehicles, appliances, and other household devices. This latter definition (def. 2) incorporates the distribution and transmission networks that are major agents in the Smart Grid. Moreover, it includes many different actors that will participate in the Smart Grid such as electric vehicles, thermostats, etc. Six categories of communications have been identified by the Department of Energy (DOE) : Advanced Metering Infrastructure (AMI), Demand Response or Demand Side Management (DR or DSM), Wide Area Situational Awareness (WASA), Distributed Energy Resources and Storage (DERS), Electric Transportation or Electric Vehicles (ET or EV) and Distribution Grid Management (DGM) [START_REF] Of | Communications Requirements of Smart Grid Technologies[END_REF]. Management of such a complex network is of great importance and this is the purpose of this paper. The electrical grid will become a smart system. Introducing smartness is all about increasing communications between local entities in order to coordinate all the elements of the global system. Indeed, it is mandatory to refer to communications when we talk about the Smart Grid. The data volume exchanged by the different actors is such that it would be impossible to transport the data without a reliable and distributed communication system. Cognitive radio [START_REF] Mitola | Cognitive radio: making software radios more personal[END_REF] [START_REF] Haykin | Cognitive radio: brain-empowered wireless communications[END_REF] could be an efficient way to transport the data without overloading the frequency spectrum [START_REF] Hossain | Smart Grid Communications and Networking[END_REF].

It has been envisaged to use different frequency bands for cognitive communications for the Smart grid [START_REF] Khan | Cognitive radio for smart grids: Survey of architectures, spectrum sensing mechanisms, and net-working protocols[END_REF], for example, free TV bands [START_REF] Qiu | Cognitive radio network for the smart grid: Experimental system architecture, control algorithms, security, and microgrid testbed[END_REF][10], in [START_REF] Yu | Cognitive radio based hierarchical communications infrastructure for smart grid[END_REF], authors consider the use of cognitive communications over unlicensed band for home communications and over licensed band for communications over longer distances. In [START_REF] Shah | A cross-layer qos-aware communication framework in cognitive radio sensor networks for smart grid applications[END_REF] the smart grid is viewed as a huge Wireless Sensors Network (WSN) in which all the elements of the smart grid com- We believe that the current centralized management of power system, is not suitable anymore for several reasons such as the increasing consumption, the introduction of distributed renewable energy (distributed geographically and in time), the customers involvement, the huge amount of data to handle, etc. We believe that, in the case of the management of such a big and complex system as the SG, the management should be distributed among all actors of this system and HDCRAM could be a management solution.

This is what we would like to highlight in this paper. Some hierarchical and distributed architectures have already been proposed for the Smart Grid. In [START_REF] Verschae | A distributed hierarchical architecture for community-based power balancing[END_REF], coordinators receive information from consumer or from the lower level coordinators and broadcast to all agent the same data which are used to make decisions. In [START_REF] Li | Distributed smart-home decision-making in a hierarchical interactive smart grid architecture[END_REF], the architecture has three levels and decision can be done at the microgid level, at the feeder level or at the substation level. In these articles hierarchical architectures are only used for demand-side management and for management pricing. In [START_REF] Dou | Multi-agent based hierarchical hybrid control for smart microgrid[END_REF], a detailed hierarchical and distributed architecture is proposed for microgrid management but this architecture hasn't been applied over a larger part of the grid as we do in this paper.

In summary, in this paper :

• A hierarchical architecture (coming from Cognitive Radio domain) is proposed for the distribution network and is used to manage consumption, intermittent production, storage and grid monitoring;

• The hierarchical position of each element of the grid (e.g. Smart Meters, Phasor Measurement Units (PMUs) or EV aggregator) is clearly identified in this paper.

• This architecture is finally implemented for Smart Home management . This architecture can also be mapped on any of the sub parts of the SG. Moreover, HDCRAM is composed by two independent paths, one for the cognitive management and the other for the (re)configuration management, that differentiates it from existing solutions. This paper is organized as follows: Section 2 describes the basics of HDCRAM, while Section 3 proposes the deployment of HDCRAM for the SG (illustrated by two clear schemes). In this section, the metrics used for both the cognitive information and the reconfiguration orders are also described for some elements of the networks. Section 4 gives a concrete example and its simulation through an agent-oriented implementation using JADE (Java Agent DEvelopment Framework). Conclusions are drawn in Section 5.

Hierarchical and Distributed Cognitive Radio Architecture Management (HDCRAM)

This section presents an architecture designed at first for Cognitive Radio purposes [START_REF] Godard | An executable meta-model of a hierarchical and distributed architecture management for cognitive radio equipments[END_REF]. Its name is Hierarchical and Distributed Cognitive Radio Architecture Management, in short: HDCRAM. The HDCRAM architecture is depicted in Figure 1. The cognitive cycle (see Figure 2) shows that a Cognitive Radio equipment or system has three main activities: a sensing activity, an intelligent activity (decision making) and a reconfiguration activity (adaptation). In HDCRAM, those activities are distinguished. The two other characteristics that are also well separated in HDCRAM, are the management architecture (level 1, 2 and 3) and the radio transceiver chain (at operator's level). Hence HDCRAM is all that has to be added to transform a non-intelligent legacy system to a smart system, whatever the system is: radio, grid, city, car, etc.

As it can be seen in Figure 1, this architecture comprises two managing sub-parts [START_REF] Moy | High-level design approach for the specification of cognitive radio equipments management apis[END_REF]:

• A cognitive part featuring Cognitive Radio Management Units (CRMu):
a CRMu exchanges information from the lower level to the upper level Communication can also exist from a CRMu to a ReMu at the same level. Figure 1 shows also a three-levels architecture:

• Level 1 (L1): is composed of the couple L1 CRM/L1 ReM and is the general manager of the system and is unique;

• Level 2 (L2). Level 2 is composed of several couples (L2 CRMu/L2 ReMu).
They provide L1 with an abstraction of the lower level L3 by not sending to L1 all the information they receive from L3 but only the necessary information. Then L2 acts as a compression level in order to decrease the amount of data to be handled in the management architecture;

• Level 3 (L3). Level 3 is composed of several couples (L3 CRMu/L3 ReMu).

Each of these couples is associated to an operator. L3 ReMu is the entity that is able to reconfigure its associated operator, and L3 CRMu is used to process the data that come from the operator.

• Operators: an operator is a modifiable component of the managed system (for example a filter or a signal to noise ratio sensor in the radio context).

The intelligence is distributed because the real-time constraints would be difficult to handle if all the decisions were made in a centralized way. This hierarchical and distributed architecture allows to make decisions at different level [START_REF] Moy | High-level design approach for the specification of cognitive radio equipments management apis[END_REF]. A local and simple decision is made at level L3 (little cycle in Figure 3). If the decision is more complex or involves other operators, the decision is made at level L2 (Figure 4 ). Lastly, if many operators are involved, the decision is made at level L1 (Figure 5). For example, Figure 6 zooms on the sub-part of a smart radio transceiver that could change the modulation when the Signal to Noise Ratio is below a given threshold [START_REF] Palicot | Radio Engineering: From Software Radio to Cognitive Radio[END_REF]. In this context, when L3 CRMu management unit of SNR sensor detects that the threshold is crossed, L3 CRMu of the SNR sensor can not directly affect the configuration of the mapping operator through its L3 ReMu. As a consequence, L3 CRMu of the SNR sensor has to refer to L2 CRMu of L2 Modulation management that can reache Mapping management L3 ReMu to manage the reconfiguration of the operator (medium decision cycle is used).

In order to help the reader really understand this HDCRAM architecture, Figure 7 provides an example of the OSA (Opportunistic Spectrum Access) challenge. OSA is a futuristic spectrum access scheme aiming to mitigate radio-frequency spectrum congestion. With OSA, future radio devices will be more agile in frequency and then should be able to change their transmitting frequency dynamically. It consists at giving new rights to so called Secondary Users (SUs) that could use the spectrum let vacant from time to time by licensed users, e.g. Primary Users (PUs). This scenario is only acceptable if SUs never interfere with PUs. A SU runs a spectrum holes sensor that senses its radio environment and sends an information about the presence or absence of a PU communicating in the targeted channel to its corresponding L3 CRMu. If a channel is free, e.g. not used by a PU, there is a hole in the spectrum and a SU can use it. If a PU uses the channel, the SU cannot use this channel and should seek for other channel opportunities and consequently should reconfigure its radio parameters such as carrier frequency. Here is the detail of the reconfiguration process that occurs inside the cognitive equipment from the sensor to the reconfiguration action. In this scenario, The L3 CRMu in charge of the sensor can't make the decision but should send a piece of information to its corresponding L2 CRMu.

Finally, this last entity concludes that the channel is no longer available.

In this scenario, L1 CRM has been designed to provoke a changing of the carrier frequency in order to explore another radio channel, thus L2 and L3 in charge of spectrum sensing can't make the decision. Consequently, the sensing result of "Spectrum Holes Sensor" is transmitted through L3 CRMu and then L2 CRMu to inform L1 CRM, which then sends a reconfiguration order to L1 ReM. L1 ReM sends a reconfiguration order to the L2 ReMu in charge of channel selection and to the L2 ReMu in charge of Common Pilot Channel (CPC) necessary to coordinate SUs that are communicating together, so that they can exchange their frequency parameters. Finally, the SU can check if the new band is free.

In order to have a high-level specification of this architecture, HDCRAM was modeled by using the Unified Modeling Language (UML) language. The use of an executable programming language, such as C, C++, Java, etc., in combination with UML and a metamodel 2 , simulates the architectures that are modelled [START_REF] Moy | High-level design approach for the specification of cognitive radio equipments management apis[END_REF].

The HDCRAM has some particularities that make it interesting for Smart Grid management, the main ones are :

• Two totally independent paths for the cognitive management and for the (re)configuration management.

• HDCRAM is a model of an architecture independent of the process-2 A metamodel is used to make the link between a diagram and a code template. • "HDCRAM model" is implemented in a specific "HDCRAM implementation" for a given smart scenario.

• Three levels of possible decisions, which corresponds to three different cognitive cycles in the architecture.

• HDCRAM model could be applied to any complex systems. The rules and the connections between the "boxes" of the architecture help the designer of the smart system to deploy smart scenarios and to specify the different elements necessary to implement this scenario. Therefore it could be applied to the entire SG comprising energy production (monitoring renewable energy,etc.) and not limited to home environment.

• An "HDCRAM implementation" could be emulated, simulated in order to predict the behavior of a system following a given smart scenario.

• An "HDCRAM implementation" could embed any decision making algorithms.

Application of Hierarchical and Distributed Cognitive Architecture for the Smart Grid

HDCRAM has been proposed within the framework of Cognitive Radio but it could be applied to any big and complex system with self-adaptation requirements, such as the Smart Grid. In brief, HDCRAM is a management tool used at first to represent a Cognitive Radio system but can be generalized to represent and manage complex and heterogeneous self-adaptive systems. HDCRAM could be an interesting solution to manage efficiently the grid. Hierarchical and distributed architectures already exist for the power system, but it is often related to the transmission part of the grid. With this HDCRAM approach, these hierarchical and distributed architectures can be generalized to the entire grid or to any sub part. The distribution part of the grid is the one whose management needs will increase the most in the coming years, that's why we maps the HDCRAM on this sub part. This gives moreover the opportunity to merge a lot of the Smart Grid aspects into a unique coherent system.

Mapping on the distribution network

Figure 8 and 9shows two sides of the same representation, cut for clarity purposes, of HDCRAM deployment for the Smart Grid. For many re-searchers and industrial stakeholders of the Smart Grid community, such as [START_REF] Moslehi | A reliability perspective of the smart grid[END_REF],improvements must be performed in the distribution network to enable the Smart Grid. Indeed, a lot of topics currently studied are related to the distribution network: demand side management programs that are enabled by distribution companies, the integration of renewable energy and storage facilities. This is, most of the time, located at the distribution level, and dealt with distribution networks monitored by PMUs, etc. This is the reason why the distribution data center is the level one (L1) in Figure 8.

Level two (L2) divides the distribution management into its three main challenges:

• The network monitoring, being represented by phasor data concentrator (fig. 8) (PDC);

• The consumers Participation (DSM vision), being represented by the home smart meter and the industrial smart meter (fig. 8).

• The integration of renewable energy and storage facilities, which is represented by the wind farms' monitoring, the PV farms' monitoring and the EV aggregator (fig. 9).

It is important to notice that, in those Figures, the energy grid and the management architecture are well separated.

With this proposed architecture, distribution companies are relieved from a part of the management. Indeed, local problems are solved locally thanks to the distributed management. If, for instance, the energy consumption of customers, in a certain area, is too large and there is therefore a pressure on the grid (information received from a PDC), smart meters will first try to solve the problem by adapting or shifting some loads or propose customers to do so by offering them a lower electricity price (Real Time Pricing program). If the engendered load modification is not sufficient, the PDC still sends information about critical states for some lines, the distribution data center will send reconfiguration orders to EV aggregators and to renewable energy sources in order to increase electricity generation. For instance, if some wind turbines/PV panels are switched off in a wind/PV farm and there is a lot of wind/sunshine, a switching ON order will be sent.

EV aggregators would also propose interesting prices in order to encourage people to discharge their car batteries into the grid.

Description of level 2

3.2.1. Network monitoring (L2 PDC of figure 8)

A phasor data concentrator (PDC) is used to collect the data that come from all the PMUs [START_REF] Tholomier | Phasor measurement units: Functionality and applications[END_REF]. A Level three (L3) related to a PDC is necessary for a PMU to be set up on a line. In runtime operations, a PMU measures voltage, current, etc. of its dedicated line and sends information to the corresponding PDC. In order not to overflow the distribution data center and provide it with level three abstraction, the PDC can only send information to it when an issue occurs or when the situation is critical. The reconfiguration orders from level one to level two and from level two to level three are the data transmission rates of a PDC and that of a PMU respectively. Electric vehicles (EV) can offer an efficient way of electricity storage.

The possibility of offering ancillary services by using EV has already been proposed by means of using cars' batteries [START_REF] Hossain | Smart Grid Communications and Networking[END_REF]. An aggregator such as a This aggregator then proposes, via a cellular network (cars are equipped with cellular antennas), to cars around, to sell their electricity at a certain price per minute/hour. The aggregator receives benefits from having distributed information and thus having contributed to ancillary services. It is important to notice that ancillary services are not in unique direction, i.e. cars can also extract energy if an overproduction is detected. This energy extraction is at lower price, so as to attract people to charge up their vehicles.

Moreover, this EV aggregator entity acts as a charging station, i.e. it has the same role as a fuel pump.

Table 3 lists examples of cognitive information and reconfiguration orders sent in the EV aggregator context.

Consumers participation (L2 Home Smart Meter of figure 8)

The smart meter is essential to enable demand side management programs. In the coming years, a typical home will be equipped with about twenty-to-thirty smart electrical devices [START_REF] Fan | The new frontier of communications research: Smart grid and smart metering[END_REF] meter. In addition to measuring consumption the smart meter will be able to operate devices. A centralized way to process the data and control the devices, i.e. all the data from all the devices sent to the distribution data center and processed by itself, seems difficult to imagine and implement practically, whereas a distributed vision would relieve the distribution data center. Since the smart meter is the link between a home and the distribution data center, this device is at level two in the HDCRAM deployment. In this distributed vision, the smart meter has the role of home energy manager, i.e. it has to manage all the smart electrical devices (level three) that consume, produce or store electricity and sensors that sense the environment.

Between level one and level two, reconfiguration orders can be load curtailments, in such cases where a customer has subscribed for a DLC (Direct Load Control) program or set electricity prices in the smart meter. In such instances, the customer may want to curtail its load or not. The cognitive information is the total energy consumption of the home.

Level three is divided into five categories: the sheddable load, the adjustable load, devices that produce electricity, Devices that can store the electricity, sensors (of any kind). In the case where the customer is an industrial company, DSM can also be applied [START_REF] Hamid | Intelligent demand side management technique for industrial consumer[END_REF]. For example, the industrial Smart Meter can reduce power consumption of non-essential devices or switch on diesel generator during peak load hours. Thus, smart industrial management can be done in the same manner as smart home management.

More details on level 3 of consumer participation are given in the following section.

3 of consumer participation

The sheddable load

The sheddable load can be reprogrammed at another time, i.e. can have its consumption shifted to off-peak hours. For instance, a washing machine can have its consumption shifted to night hours, where prices are lower.

The adjustable load

The adjustable load that cannot have its consumption totally shifted to later hours but can adapt its power level in order to decrease a bit the total electricity consumption of the house. An air-conditioner can decrease its electricity consumption by not decreasing too much the ambient temperature for example.

Paper [START_REF] Khan | Review article: A comprehensive review of the application characteristics and traffic requirements of a smart grid communications network[END_REF] even proposes to divide the load into three categories which are sheddable, adjustable and load sensitive. The sensitive load refers to the load that cannot have its consumption shifted or adapted (or just very little adapted). But this category can be merged with the adjustable load by keeping in mind that there exists devices that have zero power level adaptation.

Electricity production

The electricity generation produced by photovoltaic panels (PV) is provided in an intermittent way and depends on the environment. A local way to control those panels is desirable, i.e. switching ON/OFF the panels depending on the user requirements or distribution lines' states.

Storage facilities

With the advent of electric vehicles (EV), the possibility of storing energy at home has emerged. If the generation of the home PV panels is too large compared to the power consumption of all the other devices, the excessive energy will be injected into the grid, which is sometimes not desirable for the sake of the distribution network's performance. This energy could instead be stored in the user car's batteries in order to either use the car or use this stored energy when electricity prices are high.

Sensors

Sensors can be used to anticipate power consumption and production.

They can be of many types irradiance sensors, wind sensors, temperature sensors, etc. For example, temperature sensors can inform the smart meter of the outside temperature. The smart meter will then decide at what time radiators must be switched ON/OFF, depending on the electricity price, grid pressure, thermal inertia of the building, etc. Wind sensors can also be useful: sensors sense the velocity of the wind, send this information to the smart meter which can decide to send it to the distribution data center and the latter can predict the electricity generation of a wind farm close to a specific location.

All the devices listed above will have to communicate with the smart meter. Wireless communications is preferred as it is less intrusive and can be installed without any alteration work. But, if the number of devices is large and if some devices produce too many electromagnetic waves (such as a microwave), some frequency bands could be congested. A proposed solution is to use devices equipped with radio reconfigurable abilities and a decision making process such as OSA, so as to select the right band and channel that can be used by each device.

A Multi-Agent Based Implementation of a Home Energy Manager

This section presents a multi-agent based implementation of a home energy manager, which can be useful to simulate the behavior of a house with respect to electricity price changes. This implementation is based on the level two smart meter and the corresponding level three parts of the architecture presented in section 3 : L3 manager of PV panels, washing machines, air conditioning, etc. The proposed implementation is restricted on these parts here, in order to have a quick simulation on a concrete case. The purpose of this simulation aims to show that the application of HDCRAM architecture provides benefits compared to the case where this architecture is not applied (without management).

Java Agent DEvelopment Framework (JADE)

In order to implement a home energy manager that includes independent entities that can communicate with each other, while maintaining a certain hierarchy (an entity cannot communicate with any other entities), the JADE platform was used [START_REF] Garcia | An intelligent agent-based distributed architecture for smart-grid integrated network management[END_REF].

JADE belongs to the agent-oriented programming (AOP) class, which is used to simulate artificial distributed systems. AOP is used to model an application as a set of autonomic, active and communicating entities called agents [START_REF] Bellifemine | Developing Multi-Agent Systems with JADE[END_REF]. With these abilities, agents can independently make decisions and send any relevant information that they computed to other agents. Since JADE is used in the distributed system, agents can be easily added without any changes in the already implemented agents (they act independently).

Simulations scenario

A multi-agent based home energy manager is implemented. The concerned sub part of figure 8 and 9 in this scenario, is the level two home smart meter and the level three entities and operators below it. The thirteen electrical appliances considered in the scenario simulated can be classified in the four categories already mentionned :

• Sheddable load : dish washer, vacuum cleaner, oven, iron, washing machine and air dryer;

• Electricity production : PV agent

• Storage facilities : two electrical vehicle;

• Adjustable load : fridge, air conditioner, ventilator and radiator;

A simplified scheme of the architecture applied to the simulation scenario is depicted in Figure 10. The aim of this simulation is to compare the instantaneous power consumed with and without the hierarchical manager. The implementation is time-dependent and models a whole day divided in 24 periods (time step = one hour, this could be changed depending on the required precision without loss of generality). At each hour h, in order to reduce the pressure on the power grid, the aim of the manager is to minimize the magnitude of the total power consumption |P tot (h)|.

The expression of the total power consumption P tot (h) may be decomposed into the sum of the power consumption of each device.

For each hour, the total power consumption can be written:

P tot (h) = P cons (h) -P prod (h). (1)
The expression of the total power consumed at hour h, P cons (h), is :

P cons (h) = Na i=1 P i adj (h) + Ns i P i shed (h) + P st (h), (2)
where N a and N s denote respectively the number of adjustable and sheddable devices, P i adj (h) the power consumption of one of the adjustable devices during hour h, P i she (h) the power consumption of one of the sheddable devices during hour h and P st (h) is the power stored into the batteries of the two electrical vehicles.This power is produced by PV panels and stored because it isn't instantaneously consumed.

The total power produced P prod (h) :

P prod (h) = P P V (h) + P del (h), (3)
where P P V (h) is the power produced by PV panels and P del (h) is the power delivered by the batteries of the two electrical vehicles where PV current is stored when it is not consumed.

Finally, if we denote ∆S(h) = P st (h) -P del (h) the power variation induced by the batteries, the total power consumption becomes :

P tot (h) = Na i=1 P i adj (h) + Ns i=1 P i shed (h) -P P V (h) + ∆S(h). (4)
It is important to note that, each of the four terms of this sum relates to one level 3 unit of the architecture.

In our scenario, each of the two electrical vehicles can store 15 kWh and the power delivered or stored by each of them is between -3000 and 3000 W.

The power produced by PV panels depends on the weather, for simulation we consider that at each hour the power produced P P V (h) is uniformly distributed between 0 W and a maximum value which depends on the time,

P P V (h) ∼ U (0, M P V (h)) .
The values of M P V (h) are referenced in Table 5.

Hours of the day [h] maximum Power generated [W]

h Each of the N a = 4 adjustable devices are listed in Table 6. When the total power consumption is high, the power consumed by those devices is reduced from their nominal power consumption to their minimal power consumption.The six sheddable devices are listed in no role in the scenario, either because they have no capability to sense and to be swithed remotely, or just because they have not been considered in the scenario. HDCRAM infrastructure is deployed where it is only necessary, depending on the use cases that should be supported by the management system.

Simulations results

Figure 12 shows the comparison between a classical utilization of the operators and the result of the HDCRAM implementation presented above for a period of 24 hours. Thanks to the proposed power management, the power consumption is zero (P tot (h) = 0) for 17 of the 24 hours of the day and the peak consumption is strongly reduced (figure 12. This gain is due to the management of storage facilities, sheddable and adjustable load:

• Batteries allow a better distribution of the produced energy, power provided by solar panels is stored during the day and used at night. This part of the management doesn't reduce the quality of service but is limited by the batteries size.

• The use of sheddable devices is better distributed over the day, this reduces significantly the peak production. Moreover, the constraint on hours of use (earliest and latest hour) limits inconvenience for the consumer.

• Reducing power consumption of adjustable devices reduces the quality of service. Consequently, the manager should use this possibility wisely. At each hour, the impact of power reduction is limited by the total adjustable load (1720 W in the scenario simulated). Therefore, with the HDCRAM implementation, the customer will have a lower electricity bill (much less energy drawn from the grid) and the distribution network will have less pressure on its grid (much less energy injected to the grid and much less variation).

Conclusion

HDCRAM is proposed as an efficient tool to model and represent the communication interactions in the Smart Grid. Moreover, this tool could be used to distribute the tasks amongst the components of the SG, and HDCRAM will link all of them together. HDCRAM can help specifying the entire management system and its framework can help with generating the simulation code first, then the implementation code. The simulation results are very promising. In future work, retailers will be include in the architecture and communication access and protocol will be studied.

Acknowledgment

Part of this work is supported by the project SOGREEN (Smart pOwer Grid for Energy Efficient small cell Networks), which is funded by the French

  municate using cognitive radio. These papers are applying Cognitive Radio to the Smart Grid from the SG communication point of view. In this paper, the application of CR is done through the application of a Cognitive Radio Manager to the Smart Grid. The Smart Grid management architecture affects the complexity of the management of cognitive radio equipments, as a consequence, the design of the architecture must take into account the communication system. In order to fulfil the communication requirements of the SG, we propose to use a Cognitive Radio Management Architecture that was initially proposed during the first CROWNCOM conference in 2006 [13]. This architecture, named HDCRAM for Hierarchical and Distributed Cognitive Radio Architecture Management, was successfully demonstrated in several Cognitive Radio (CR) scenarios such as the Dynamic Spectrum Access (DSA) [14] [15] 1 .

Figure 1 :

 1 Figure 1: Hierarchical and Distributed Cognitive Radio Architecture Management [20] over a processing chain made of reconfigurable operators. The processing chain targeted by HDCRAM was initially made of radio processing elements and we propose in this paper to target the elements of the grid.

Figure 2 :

 2 Figure 2: Simplified cognitive cycle

Figure 3 :Figure 4 :Figure 5 :Figure 6 :

 3456 Figure 3: A little cycle (includes one level L3)

Figure 7 :

 7 Figure 7: Example of HDCRAM deployment for OSA purpose

  telecommunication operator (L2) receives information about the energy re-Levels Cognitive Information ⇑ Reconfiguration Orders ⇓

Figure 10 :

 10 Figure 10: Simplified diagram of the architecture applied to the simulation scenario

Figure 11 :

 11 Figure 11: Communication paths between L3 and the L2 manager

Figure 12 :

 12 Figure 12: Simulation results of the power consumption in W in function of time, as described in section 4.3 (a) when a classical utilization of the electrical devices is done and (b) when HDCRAM is used as a management facility of the home power consumption

Table 1

 1

	lists examples of the upward (cognitive information) and down-
	ward (reconfiguration orders) data flow sent in this context.

Table 1 :

 1 Interaction between the distribution data center and a PMU3.2.2. Integration of renewable energy and storage facilitiesThis section describes the L2 of Wind farm monitoring, PV farm moni-

toring and EV aggregator of figure

9

. Renewable energy sources are difficult to handle because of the intermittent generation. The energy produced by a wind farm or a PV farm does not necessary coincide with the energy demand. The production of these renewable energy sources at time t is sometimes much larger than the consumption at this same time, which can cause a large pressure on the grid. This is why a local control of a wind farm or a PV farm is required. Indeed, if a wind/PV farm exceeds a specific production threshold (sent by the distribution data center), a monitoring system of the wind/PV farm (L2) can send orders to one or several wind turbines/PV panels (L3) to switch OFF in order to reduce the total electricity generation. Table 2 lists examples of cognitive information and reconfiguration orders sent in this context.

Table 2 :

 2 Interaction from the distribution data center to a wind turbine/PV panel quired to balance demand and production from the distribution data center.

Table 3 :

 3 Interaction from the distribution data center to cars batteries

Table 4

 4

		lists examples of cognitive information and reconfiguration or-
	ders between the level two smart meter and all the levels three related to
	this level two, as described above.
	Level 2-3	Cognitive Information ⇑	Reconfiguration Orders ⇓
		Electricity consumption	Shift consumption to time
	Sheddable	Time at which the device	t (switch OFF until time t)
		will start consuming	Switch ON/OFF the device
		Electricity consumption	Switch the device
	Adjustable Min power allowed at time t consumption to its minimum/
		Max power allowed at time t	maximum consumption
		Batteries level of the car	Change the batteries state:
	Storage	Batteries state: charge,	charge, unplug from grid, etc.
		discharged, etc.	Discharge in the grid.
	Renewable	Production level	Switch OFF/ON
	energy	of the PV panels	the PV panels
	External	Send sensing data:	Not really any
	Sensors	wind velocity, sunshine,	reconfiguration orders
		temperature, etc.	but switch OFF/ON a sensor

Table 4 :

 4 Interaction between level two smart meter and level three devices at home

Table 5 :

 5 PV agent used for the simulation

	∈ [21; 4]	30
	h ∈ [5; 7] ∪ [18; 20]	4000
	h ∈ [8; 17]	12000

Table 7

 7

	, they must be used

Table 6 :

 6 Adjustable devices agents used for the simulation

	Each hour, all L3 CRMu in charge of power production, storage and
	adjustable devices send cognitive information to the L2 CRMu (respectively
	the power produced, state of batteries and adjustable load requirements).
	This entity computes the power consumption of the hour and then sends
	configuration order (adjustable devices to start, power level of sheddable
	devices and power to deliver or store for batteries) to the L2 ReMu which
	transmits it to the corresponding L3 ReMus. Configuration orders must take
	into account the constraint of the devices pre-set in the manager. Figure
	11 shows communication paths between the L2 manager and the L3 units
	empoyed for simulations.
	Note that HDCRAM is able to support that some devices in white have

Table 7 :

 7 Sheddable devices agents used for the simulation

Table 8

 8 compares the results obtained in automatic and manual mode.It is clear that the variability and the maximum/minimum power levels in manual mode are much greater than in automatic mode.

		Peak	Peak
	Mode	consumption generation
		[W ]	[W ]
	HDCRAM	850	800
	Classical	5500	7000

Table 8 :

 8 Comparison between HDCRAM and a classical utilization

A demonstration is also available on the internet : https://youtu.be/0-MIpPwgG4E


Siamak Niroomandi

Icíar Alfaro

David González

Elías Cueto

Francisco Chinesta

Model order reduction in hyperelasticity: a proper generalized decomposition approach

Keywords: proper generalized decomposition, separated representations, asymptotic numerical methods, model order reduction, real time simulation, hyperelasticity

This paper deals with the extension of proper generalized decomposition methods to non-linear problems, in particular, to hyperelasticity. Among the different approaches that can be considered for the linearization of the doubly weak form of the problem, we have implemented a new one that uses asymptotic numerical methods in conjunction with proper generalized decomposition to avoid complex consistent linearization schemes necessary in Newton strategies. This approach results in an approximation of the problem solution in the form of a series expansion. Each term of the series is expressed as a finite sum of separated functions. The advantage of this approach is the presence of only one tangent operator, identical for every term in the series. The resulting approach has proved to render very accurate results that can be stored in the form of a meta-model in a very compact format. This opens the possibility to use these results in real-time, reaching kHz feedback rates, or to be used in deployed, handheld devices such as smartphones and tablets.

INTRODUCTION

Proper orthogonal decomposition (POD) methods have had a tremendous importance in many branches of applied sciences and engineering, where they are known under a variety of names. Starting from the pioneer works by Karhunen [START_REF] Karhunen | Über lineare methoden in der wahrscheinlichkeitsrechnung[END_REF] and Loève [START_REF] Loève | Probability Theory[END_REF] (also Lorenz [START_REF] Lorenz | Empirical Orthogonal Functions and Statistical Weather Prediction[END_REF]), POD methods, also known as principal component analysis and subspace tracking, among other names, have been applied to a wide range of problems, from chemical engineering [START_REF] Park | The use of the Karhunen-Loève decomposition for the modeling of distributed parameter systems[END_REF] to real time visualization, possibly with haptic response [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF], structural mechanics [START_REF] Idelsohn | A reduction method for nonlinear structural dynamics analysis[END_REF][START_REF] Krysl | Dimensional model reduction in non-linear finite element dynamics of solids and structures[END_REF], or virtual surgery [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF][START_REF] Taylor | A reduced order explicit dynamic finite element algorithm for surgical simulation[END_REF][START_REF] Taylor | A reduced order finite element algorithm for surgical simulation[END_REF].

Proper orthogonal decomposition rely upon the computation of the solution to surrogate problems, the so-called snapshots. The most typical structure of these snapshots, obtained as an eigenvalue problem through a minimization process, is then used as a basis for subsequent problem solving. Among the problems related to this approach, one can cite the need for an important number of snapshots. These snapshots should be generated by solving problems somewhat similar to the one at hand, but maybe for slight changes in geometry or boundary conditions. The so-called a priori POD [START_REF] Ryckelynck | A priori model reduction method for the optimization of complex problems[END_REF][START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF][START_REF] Ryckelynck | On the a priori model reduction: overview and recent developments[END_REF] constitutes one of the first attempts to avoid the costly generation of snapshots by starting from scratch and adding Krylov subspaces to the basis of the reduced model.

Another difficulty on the use of POD-based model order reduction is related to the need of recomputing the original system of equations in non-linear problems (even if it will not be solved) to update tangent stiffness matrices. This prevents the method to provide as competitive costs as one would expect from the limited size of the resulting systems of equations. This problem has generated recently quite an active area of research. As notable contributions, one can cite the empirical interpolation method [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF] or its discrete counterpart [START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF]. A somewhat related approach can be found in [START_REF] Nguyen | A 'best points' interpolation method for efficient approximation of parametrized functions[END_REF].

On the other hand, non-linear POD methods (those designed for data lying on a manifold) are always difficult to develop. See [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF][START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF][START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF][START_REF] Tenenbaum | A global framework for nonlinear dimensionality reduction[END_REF] to name a few. Interpolation among reduced models continues to be an issue too [START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF][START_REF] Amsallem | An interpolation method for adapting reduced-order models and application to aeroelasticity[END_REF].

Proper generalized decomposition (PGD) methods, on the contrary, arose recently as a generalization of POD techniques. PGD roots can be traced back to the pioneer work by P. Ladeveze on the LATIN method [START_REF] Ladeveze | Nonlinear Computational Structural Mechanics[END_REF] and, particularly, the so-called radial loading approximation scheme, in which a separated space-time approximation of the displacement in structural mechanics problems was used. Independently, the method was re-invented in the framework of problems defined in high dimensional state spaces [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: transient simulation using space-time separated representations[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF]. It was then soon realized that PGD methods can be considered as a generalization of POD, in which the basis are computed on the fly without no previous snapshot. Instead, the essential field is approximated as a finite sum of separable functions, very much similar to the space-time structure of the radial approximation within the LATIN method. For recent surveys on PGD methods, the reader is referred to [START_REF] Chinesta | Recent advances in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF][START_REF] Ladeveze | The Latin multiscale computational method and the proper generalized decomposition[END_REF]. Some of the nice properties of POD are lost, however. For instance, optimality of the POD basis is no more guaranteed, although convergence properties have been demonstrated recently [START_REF] Bris | Results and questions on a nonlinear approximation approach for solving highdimensional partial differential equations[END_REF], and error estimators have also been proposed [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF][START_REF] Ladeveze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF].

When compared with the existing model order reduction techniques, such as the ones based on POD, for instance, PGD provides the modes a priori, that is, without the need to compute snapshots of the complete model. These snapshots would be extremely costly to generate for a problem defined in high dimensional phase spaces, where the sampling of the space is not obvious. With respect to the technique introduced in [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF], on the contrary, it is stated that 'deformation basis generation is a hard open problem in solid mechanics, and there exist no algorithms for automatic proven-quality global deformation basis generation under general forcing'. This is precisely where PGD comes into play. The method here proposed provides with a suitable set of reduced basis a priori, without the need to compute costly snapshots. In [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF], two different methods were chosen to generate the basis; one based upon standard modal analysis that includes derivatives of the eigenmodes. The second one is basically a POD procedure in which the modes are extracted from a linear simulation and successively enriched off-line. However, although in this reference, it is stated that the method is general and can be applied to constitutive models other than Kirchhoff-Saint Venant (which is of low practical interest because of its well-known instabilities in compression due to lack of convexity of the strain energy functional), we believe that it is difficult to generalize for state-of-the-art constitutive hyperelastic models. Up to our knowledge, it has not been generalized yet.

One key aspect in the rapid development of PGD methods is related to the fact that PGD can be seen as both an efficient technique for high dimensional problems and as an a priori model order reduction technique. This opened the door to re-interpreting parametric problems as high dimensional ones, just by considering parameters as new dimensions of the state space of the model [START_REF] Ch | Methodological approach to efficient modeling and optimization of thermal processes taking place in a die: application to pultrusion[END_REF][START_REF] Ch | Proper generalized decomposition-based dynamic data-driven control of thermal processes[END_REF][START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the Proper generalized decomposition[END_REF].

This new point of view offered new and interesting insights in many branches of applied science and engineering. It is now envisaged how by looking at a parametric problem from a multidimensional framework gives rise to important savings in simulation time. Furthermore, it has been observed how important advantages could be obtained by formulating standard problems as high dimensional ones. For instance, formulating plate and shell problems as 3D ones but with 2D computational complexity [START_REF] Bognet | Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity[END_REF] has revealed very interesting insights on a general finite element (FE) formulation beyond Kirchhoff and Reissner-Mindlin plates.

In this work, we explore the possibilities of PGD methods applied to fast (real time) simulation of hyperelastic solids. For a related reference, the interested reader can also consult [START_REF] Relun | A model reduction technique based on the PGD for elastic-viscoplastic computational analysis[END_REF]. It will be shown how a multidimensional formulation of the problem, in which the displacement is expressed as to be dependent on both the physical coordinates and the position and orientation of the applied loads (thus, defined in R 9 ), opens the door to simulations with feedback rates on the order of 1 kHz. The developed strategy comprises an off-line computation strategy, in which a highdimensional problem is solved. This solution provides in fact a sort of metamodel that can be stored in a very compact form. Then, an online simulation strategy is developed that solves the metamodel at impressive feedback rates, even on handheld devices.

But the price to pay in this case is the complex non-linear structure of the weak form of the problem, as will be shown in Section 2. Consistent linearization of PGD weak forms continues to be an issue in our community. Previous works include explicit linearizations [START_REF] Niroomandi | Real time simulation of biological soft tissues: a PGD approach[END_REF]. In this work, we introduce a new approach by combining a PGD formulation of the problem with an asymptotic expansion of the displacement. The use of asymptotic expansions in the search of complex nonlinear equilibrium paths of structures is due to Potier-Ferry and coworkers [START_REF] Abichou | Asymptotic numerical method for problems coupling several nonlinearities[END_REF][START_REF] Cao | An improved iterative method for large strain viscoplastic problems[END_REF][START_REF] Cochelin | Asymptotic-numerical methods and Padé approximants for non-linear elastic structures[END_REF][START_REF] Cochelin | The asymptotic numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF]. This gives rise to a formulation in which only one tangent operator is present, the same for all the terms of the expansion, and where the non-linearities are translated to the right-hand-side (force) term of the matrix form of the problem. The resulting technique provides results, as will be shown, at kHz feedback rates for general hyperelastic models. This constitutes, in our opinion, a major advancement in the state of the art of reduced order models. Up to our knowledge, only techniques working with Kirchhoff-Saint Venant model were able to obtain such performance rates.

In Section 3, we review the basics of asymptotic numerical methods (ANM) in their standard form, whereas in Section 4, we develop the combined PGD-ANM formalism and study how it provides a very convenient form for non-linear path following of hyperelastic solids and structures. In Section 5, we study the performance of the proposed technique by employing some classical benchmark tests and show how the resulting formalism can be employed advantageously to perform real-time simulation.

FORMULATION OF THE PROBLEM UNDER THE PGD POINT OF VIEW

As mentioned before, one of the main advantages of using PGD for model order reduction purposes relies in the possibility of rewriting the model in a multidimensional framework. Consider, for simplicity, the static equilibrium equations of a general solid under small strain assumptions,

r C b D 0 in , ( 1 )
where b represents the volumetric forces applied to the body, subjected to the following boundary conditions:

u D N u on u (2) n D N t on t . ( 3 )
The standard weak form of the problem is obtained after multiplying both sides of Equation ( 1) by an admissible variation of the displacement, u , and integrating over the domain . To fully exploit the characteristics of PGD methods, following [START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the Proper generalized decomposition[END_REF], we convert the equilibrium problem given by Equations ( 1)-( 3) into a parametric one, by considering, again for the sake of simplicity in the exposition, that the load N t is punctual, unitary, vertical, and is applied at a position s 2 t that acts as a parameter in the formulation. Note that this alternative formulation of the problem gives rise to a new one defined in R 6 , because u D u.x, s/ 2 N , where N Â t represents the portion of the boundary where the load can be applied.

If inertia terms in Equation (1) were non-negligible, a possible approach that considers initial conditions as additional parameters has also been explored in [START_REF] Gonzalez | Proper generalized decomposition based dynamic data driven inverse identification[END_REF].

In this spirit, an alternative (doubly) weak form of problem (1)-(3) consists in finding the displacement u 2 H 1 . / L 2 . N / such that for all u 2

H 1 0 . / L 2 . N / Siamak5, Z N Z .r s u / T d d N D Z N Z t2 .u / T tdd N , ( 4 )
where D u [ t represents the boundary of the solid, divided into essential and natural regions, and where t D t1 [ t2 , that is, regions of homogeneous and non-homogeneous, respectively, natural boundary conditions. In turn, we assume t D e ´ı.x s/, where ı represents the Dirac delta function and e ´represents the unit vector along the ´direction, in this case. Note that a general-form load term does not include further complexity into this formulation.This Dirac delta term should be regularized for computation purposes and approximated by

t j n X iD1 f i j .x/g i j .s/ (5)
by simply performing a singular value decomposition of the load, for instance. As explained before, PGD methods assume a separated representation of the unknown field (here, the displacement). This is the key ingredient of the method that allows solving efficiently high dimensional problems. The so-called curse of dimensionality associated to mesh-based solution of high dimensional problems is thus avoided by solving a sequence of low-dimensional problems in separated form. This was in fact the key ingredient in the space-time separated representation of the displacement in the LATIN method [START_REF] Ladeveze | Nonlinear Computational Structural Mechanics[END_REF]. The PGD approach to the problem is characterized by the construction, in an iterative way, of an approximation to the solution in the form of a finite sum of separable functions. Assume that we have converged to a solution, at iteration n, of this procedure,

u n j .x, s/ D n X kD1 X k j .x/ Y k j .s/, ( 6 )
where the term u j refers to the j th component of the displacement vector, j D 1, 2, 3.

The following term of this approximation, the (n C 1)th one, will be similar to

u nC1 j .x, s/ D u n j .x, s/ C R j .x/ S j .s/, ( 7 )
where R.x/ and S .s/ are the sought functions that improve the approximation. In this same way, the admissible variation of the displacement will be given by

u j .x, s/ D R j .x/ S j .s/ C R j .x/ S j .s/. ( 8 )
At this point, several options are at hand so as to determine the new pair of functions R j and S j . The most frequently used, due to both its ease of implementation and good convergence properties, in general, is a fixed-point alternating directions algorithm in which functions R j and S j are sought iteratively. We describe briefly the implementation of this algorithm.

Computation of S.s/ assuming R.x/ is known

In this case, following standard assumptions of variational calculus, we have

u j .x, s/ D R j .x/ S j .s/, ( 9 )
or, equivalently, u .x, s/ D R ı S , where the symbol 'ı' denotes the so-called entry-wise Hadamard or Schur multiplication for vectors. Once substituted into Equation ( 4), gives

Z N Z r s .R ı S / W C W r s n X kD1 X k ı Y k C R ı S ! d d N D Z N Z t2 .R ı S / m X kD1 f k ı g k ! dd N , (10)
or equivalently (we omit obvious functional dependencies)

Z N Z r s .R ı S / W C W r s .R ı S /d d N D Z N Z t2 .R ı S / m X kD1 f k ı g k ! dd N Z N Z r s R ı S R n d d N , (11)
where

R n represents R n D C W r s u n . ( 12
)
All the terms depending on x are known, and hence we can compute all integrals over and t2 (support of the regularization of the initially punctual load) to derive an equation to compute S .s/.

Computation of R.x/ assuming S.s/ is known

Equivalently, in this case, we have

u j .x, s/ D R j .x/ S j .s/, (13)
which, once substituted into Equation ( 4), gives

Z N Z r s .R ı S / W C W r s n X kD1 X k ı Y k C R ı S ! d d N D Z N Z t2 .R ı S / m X kD1 f k ı g k ! dd N . (14)
In this case, all the terms depending on s (load position) can be integrated over N , leading to a generalized elastic problem to compute function R.x/.

This simple algorithm renders generally excellent convergence properties (see [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF] and references therein).

While fixed-point approaches (see [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF] and references therein) render generally excellent results, it is generally very difficult to perform consistent linearizations on Equations ( 11) and [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF]. In [START_REF] Niroomandi | Real time simulation of biological soft tissues: a PGD approach[END_REF], some explicit algorithms have been tested to overcome this problem that generally work well. By this, we mean that they provide errors, for reasonable pseudo-time step sizes, of the same order of those obtained by existing POD-based approaches such as [START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF][START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF]. This level of error is generally higher than that of common engineering practice but is nevertheless acceptable in the real-time simulation community [START_REF] Bro-Nielsen | Real-time volumetric deformable models for surgery simulation using finite elements and condensation[END_REF].

In what follows, a method that combines asymptotic expansions of the displacement and the second Piola-Kirchhoff stress tensor to avoid complex stiffness matrix update procedures is reviewed. The method, due to Potier-Ferry and coworkers [START_REF] Abichou | Asymptotic numerical method for problems coupling several nonlinearities[END_REF][START_REF] Cao | An improved iterative method for large strain viscoplastic problems[END_REF][START_REF] Cochelin | Asymptotic-numerical methods and Padé approximants for non-linear elastic structures[END_REF][START_REF] Cochelin | The asymptotic numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF][START_REF] Yvonnet | A model reduction method for the post-buckling analysis of cellular microstructures[END_REF], has been used in a variety of structural mechanics problems to follow complex equilibrium paths of solids and structures.

A BRIEF REVIEW OF THE ANM

For completeness and following closely [START_REF] Yvonnet | A model reduction method for the post-buckling analysis of cellular microstructures[END_REF], we review here the basics of the ANM applied to the simplest non-linear case, namely, Kirchhoff-Saint Venant solids (i.e., linear elastic solids undergoing large strains). Under a Lagrangian frame of reference, we consider the displacement as given by

x D X C u. ( 15
)
Following the notation in [START_REF] Cochelin | The asymptotic numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF], we consider a linear and a non-linear term for the Green-Lagrange strain tensor, E , in the form

E D 1 2 F T F 1 Á D l .u/ C nl .u, u/, (16)
where F D r u C I is the gradient of deformation tensor and

l .u/ D 1 2 r u T C r .u/ , ( 17a
) nl .u, u/ D 1 2 r u T r .u/. (17b)
Hyperelastic materials are based on the assumption of a particular strain-energy function, ‰. Then the second Piola-Kirchhoff stress tensor S can thus be obtained by

S D @‰ @E , (18)
that is a symmetric tensor and is related to the first Piola-Kirchhoff stress tensor, P, by P D F S . The equilibrium equation stated in the reference configuration is

r P C B D 0 in 0 , ( 19
)
in which B is the body force. The boundary conditions of the body are defined by (we do not consider the case of follower loads for simplicity)

u.X / D N u on u , PN D t on t , ( 20
)
where N is the unit vector normal to D @ 0 , N t is an applied traction, and is a loading parameter, equivalent to a pseudo-time and ranging from 0 to 1. The weak form of the problem is then given by

Z 0 S W E d D Z t N t u d 8u 2 H 1 . /, (21)
where in the aforementioned equation, E is expressed by

E D 1 2 h F T r .u / C r .u / T F i D l .u / C nl S .u, u /, (22)
where, in turn, nl S .u, u / is defined by

nl S .u, u / D nl .u, u / C nl .u , u/. (23)
The Kirchhoff-Saint Venant model is characterized by the energy function given by

‰ D 2 .tr.E // 2 C E W E , ( 24
)
where and are Lame's constants. The second Piola-Kirchhoff stress tensor can be obtained by

S D @‰.E / @E D C W E (25)
in which C is the fourth order constitutive (elastic) tensor.

Although it is well known that the Kirchhoff-Saint Venant model is unstable under compression, and thus of limited importance in engineering applications, we have considered it here for the sake of simplicity in the formulation and because it is still in use in fields such as real-time simulation of surgery and rendering in general, where it is among the state of the art models [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF][START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF][START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF].

Under these assumptions, the ANM is based upon expanding the displacement in the neighborhood of each material point in terms of a control parameter 'a'. This expansion is developed in the neighborhood of a known equilibrium point u m I S m I m Á at step m, and the series is truncated at order N . To simplify the resulting expressions, also the second Piola-Kirchhoff stress tensor and the load parameter are expanded in series prior to their introduction in the weak form of the problem,

where u p , S p , p Á are unknowns. Previously, u mC1 .a/, S mC1 .a/, mC1 .a/ Á represents the solution along a portion of the loading curve. Noteworthy, the behavior of the solid is described continuously with respect to 'a'. The introduction of Equation [START_REF] Ladeveze | The Latin multiscale computational method and the proper generalized decomposition[END_REF] into Equations ( 21) and ( 25) leads to a series of linear problems with the same tangent operator, thus avoiding the burden associated with stiffness matrix updating in the Newton-Raphson scheme.

In general, any variable can be expanded in terms of a, so, for instance, the series expansion of E .u/ gives

E mC1 .a/ D l .u / C nl S .u , u m / C N X pD1 a p nl S .u , u p /. (27)
As can be noticed, E mC1 includes terms in both u and u . This is not surprising, because we are computing a correction to u m . The series expansions of S gives in turn

S mC1 .a/ D C W E mC1 .a/ D C W 2 4 nl .u m , u m / C l .u m / C N X pD1 a p l .u p / C nl S .u m , u p / C p 1 X iD1 nl .u i , u p i / ! 3 5 , (28)
and at order p, we obtain

S p D C W ´ l .u p / C nl S .u m , u p / C p 1 X iD1 nl .u i , u p i / μ . ( 29
)
Introducing the asymptotic expansion into Equation ( 25) results in

Z 0 8 < : 0 @ S m C N X pD1 a p S p 1 A W 0 @ l .u / C nl S .u m , u / C N X pD1 a p nl S .u p , u / 1 A 9 = ; d D 0 @ m C N X pD1 a p p 1 A ‰ ext .u /, (30)
with ‰ ext .u / D R t t u d . Introducing Equation (29) into Equation ( 30) and identifying terms with the same power of 'a' results in a successive series of linear problems, which at order p, .p D 1, : : : , N / takes the form

L.u , u mC1 / D p ‰ ext .u / C F nl p .u / (31)
with

L.u , u mC1 / D Z ® S m W nl S .u m , u / C OE l .u / C nl S .u p , u / WCW OE l .u p / C nl S .u m , u p / ¯d (32)
and where F nl p .u / is equal to zero at order one, and at order p, it can be calculated as

F nl p .u / D Z ´p 1 X iD1 S i W nl S .u p i , u / C p 1 X iD1 OE nl .u i , u p i / WC WOE l .u / C nl S .u m , u / μ d (33)
Discretization of Equation ( 31) by using FEs leads to a sequence of linear problems in the form [START_REF] Cochelin | The asymptotic numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF] Order 1

´K t u 1 D 1 f u T 1 u 1 C 2 1 D 1, (34)
Order

p ´K t u p D p f C f nl p .u i / i < p u T p u 1 C p 1 D 0, (35)
where K t denotes the tangent stiffness matrix associated with Equation ( 32), common to the problems at different orders p. It is the same as the one applied in a classical iterative algorithm such as Newton-Raphson (in the first iteration). In the aforementioned equation, u p is the discretized form of the displacement field at order p, f is the loading vector, and f nl p represents the discretized form associated with F nl p .u / in Equation ( 33), which at order p only depends on the values of u i , i < p.

A COMBINED PGD-ANM APPROACH TO HYPERELASTICITY

The combination of the two previously introduced tools, namely the PGD approach to parametrized problems (in this case the position of the load is the parameter) and the ANM for a consistent linearization of the weak form of the problem gives rise to a particularly useful formulation. In it, the displacement field is approximated as a series expansion around the last equilibrium point, whereas each term of the series is considered to be approximated by a finite sum of separated functions,

u mC1 D u m C a n 1 X iD1 F 1 i ı G 1 i C a 2 n 2 X j D1 F 2 j ı G 2 j C : : : C a p n p X lD1 F p l ı G p l . ( 36
)
This gives rise to a series of problems of the form ( 34)-( 35), within which a traditional PGD problem is solved for functions F j i .x/ and G j i .s/ at each order of the expansion j . No further modification of the method is necessary, resulting in a series of standard PGD problems that can be solved by employing any of the available non-linear solvers. Here, as in many of our previous works, we have employed a fixed point algorithm, similar to the one sketched in Sections 2.1-2.2 before.

We study the behavior of the proposed technique by means of a series of benchmark problems in Section 5.

NUMERICAL EXAMPLES

Kirchhoff-Saint Venant material

We will consider two sets of examples. The first one is composed by three different benchmark tests over Kirchhoff-Saint Venant materials. Despite its simplicity and well-known limitations, Kirchhoff-Saint Venant approaches are very useful in the field of real-time simulation, because they provide a good compromise between realism in the deformation and computational cost [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF][START_REF] Taylor | A reduced order finite element algorithm for surgical simulation[END_REF].

Pinched cube.

A unit cube modeled by 3 3 3 nodes and a tetrahedral mesh is considered. Young's modulus of 1 MPa and Poisson's coefficient of 0.25 are assumed. The cube is loaded by a vertical force (0.01 N) acting at any point of the top face. Results obtained with the before presented PGD-ANM approach are compared with traditional FE approaches, solved by standard Newton-Raphson linearization strategies. For one particular position of the load (one corner), the load-displacement curve is shown in Figure 1, whereas its deformed configuration is shown in Figure 2.

In general, as was the case in previous references such as [START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF], expansions up to order six are enough to obtain a reasonable accuracy, despite the fact that in the ANM community, much higher order expansions are usually employed on the order of 15 [START_REF] Cochelin | Asymptotic-numerical methods and Padé approximants for non-linear elastic structures[END_REF]. The number of modes (separated functions) employed at each expansion term in this example was 17, 17, 11, 2, 1 and 1, respectively, for terms of orders one to six. The accuracy of the approach is noteworthy, despite the fact that terms 5 and 6 of the expansion (solution for the sixth order expansion is not depicted in Figure 1, because it is practically indistinguishable from that of order five) are obtained with only one couple F .x/ ı G .s/, that is, n 5 D n 6 D 1 in Equation [START_REF] Abichou | Asymptotic numerical method for problems coupling several nonlinearities[END_REF].

Beam bending.

In this case, we consider the problem of beam bending. This simple test is among the most popular ones in the field of real-time simulation, because it readily shows a great divergence from physical results if a poor formulation is used [START_REF] Allard | SOFA -an open source framework for medical simulation[END_REF]. In fact, if linear elasticity formulations are employed, a great gain of volume is observed, leading the observer to perceive a clearly non-physical result, see Figure 3.

In this case, we consider a squared cross-section beam, with 40 40 400 mm. Young's modulus was assumed to be 209, 000 MPa, whereas Poisson coefficient was set to 0.3. A load of 10 6 N can be applied at any point of the upper face of the beam. Under these conditions, a comparison was established between the tip displacement obtained for the load position at the beam rightmost side with the value obtained by using standard FE analysis and Newton-Raphson iterations to solve the non-linear equations.

Results for different expansion order are shown in Figure 4. For a fifth order expansion, almost no difference can be perceived between the FE result and the reduced PGD-ANM result. For this example, the number of separated functions necessary at each expansion order was 162, 102, 79, 161, 24, and 8, respectively. As stopping criterion, we employed a relative one, which stops the fixed point algorithm if the new functional pair contributes less than 10 4 times the initial pair of functions, and a general one in which the algorithm is stopped if the modulus of the new functional pair is less than 10 15 .

In general, results are below 5% error with respect to the target FE solution of the problem. In any case, a higher number of modes for each term, or a higher number of terms in the expansion, has been shown to converge to the right solution. It is true, in general, that the number of off-line computed modes depends on the error estimator considered (see the works by Huerta [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF] or Chamoin and Ladeveze [START_REF] Ladeveze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF]). At present, there are no efficient and robust error estimators applicable to general multiparametric non-linear models. For this reason, we computed many terms (remember that this calculation is carried out off-line and only once) and then we compared the solution for some particular choices of the parameters with the one obtained with FEs to check the convergence and solution accuracy.

5.1.3. Palpation of the liver. One of the most typical examples in real-time applications is that of liver palpation in a laparoscopic virtual surgery procedure. The liver is the biggest gland in the human body, after the skin. Liver geometry has been obtained from the SOFA project [START_REF] Allard | SOFA -an open source framework for medical simulation[END_REF] and post-processed to obtain a mesh composed by 2853 nodes and 10,519 tetrahedra, whose geometry is shown in Figure 5. The anterior surface of the liver is considered free, whereas the posterior one was assumed to be supported over different organs (it is connected to the diaphragm by the coronary ligament, for instance). The inferior vena cava travels along the posterior surface, and the liver is frequently assumed clamped a that location. Although the assumed boundary conditions are not strictly correct from a physiological point of view, our main interest is to show that the model can be solved under real-time constraints with reasonable accuracy.

Here, the human liver is considered as a Kirchhoff-Saint Venant material with E D 0.17 MPa and D 0.48 [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF]. A vertical load of 5 N is considered at any point of the anterior surface of the liver. Again, for comparison purposes, results are evaluated at a particular position of the load and compared with those obtained with a traditional FE analysis and standard Newton-Raphson iterative schemes. It can be noticed how the sixth order expansion gives very reasonable results if compared with that of the FE model (Figure 6). The vertical displacement field for that particular position of the load is depicted in Figure 7. To check the overall behavior of the technique, functions F i and G i (Equation ( 36)) for i D 1, 2, 3, and 20 (terms of first and fourth order expansions) are shown in Figures 8-11.

-0.5 -1 -1.5 -2 -2.5 -3 (b) Z X Y F 3.5 3 2.5 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5 (c) Z X Y F 9 8 7 6 5 4 3 2 1 0 -1 -2 -3 -4 -5 -6 -7 -8 (d)

Neo-Hookean behavior

Extension of the before presented technique to neo-Hookean materials [START_REF] Bonet | Nonlinear Continuum Mechanics for Finite Element Analysis[END_REF] is in principle straightforward, although lengthy (this is not important, in fact, because the parametric calculation is carried out off-line). The major difference with the Kirchhoff-Saint Venant model is the presence of material non-linearities, in addition to the geometrical ones. A POD-ANM approach has been suggested in a previous work of the authors, see [START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF], and is here extended to the before presented PGD framework.

The compressible neo-Hookean model is characterized by a strain energy function given by

‰ D 2 .tr.C / 3/ ln J C 2 .ln J / 2 , ( 37
)
where and are Lame's constants and C D I C 2E is the right Cauchy-Green strain tensor. The second Piola-Kirchhoff stress tensor can be obtained by In this case, the intricate expansion procedure becomes easier if we identify, as in [START_REF] Cao | An improved iterative method for large strain viscoplastic problems[END_REF], the asymptotic expansion with a Taylor series of the variables of interest, denoted by U .a/, in the vicinity of a D 0. Truncating at order N .

S D @‰.E / @E D .I C 1 / C .ln J /C 1 . (38)
U .a/ D U 0 C N X pD1 U p a p , (39)
where U 0 D U .0/ and In this case, as in [START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF], we have selected the following variables to perform the expansion,

U p D 1 pŠ d p U da p aD0 . ( 40
U .a/ D 0 B B B B B B B B @ u.a/ S .a/ J 2 C 1 .a/ ln p J 2 .a/ 1 J 2 .a/ .a/ 1 C C C C C C C C A . ( 41
)
By performing the substitution of the aforementioned variables into the weak form of the problem (Equation ( 21)), we arrive to a problem entirely similar to that in Equations ( 34) and [START_REF] Niroomandi | Real time simulation of biological soft tissues: a PGD approach[END_REF]. The entire details are provided, for completeness, in Appendix 1.

Neo-Hookean beam under bending.

We reproduce here the problem in Section 5.1.2 but considering a neo-Hookean constitutive model. The deformed configuration of the beam at expansion orders one to four is depicted in Figure 12. Note how the first order expansion (linear approach) shows a tremendous gain in volume that renders the simulations clearly non-physical. Again, expansions up to orders four to six were judged sufficient to obtain a good approximation to the reference solution (Figure 13). Obviously, higher accuracy can be obtained by increasing even more the expansion order. Up to p D 15 is a typical value of the expansion order in the ANM literature.

Palpation of a neo-Hookean liver.

The same procedure has been applied to the problem in Section 5.1.3 but now considering neo-Hookean behavior. The neo-Hookean law in Equation [START_REF] Cao | An improved iterative method for large strain viscoplastic problems[END_REF] has been now particularized to E D 0.17 MPa and D 0.48. As in previous examples, a PGD-ANM solution has been obtained and compared with a standard FE solution at a particular node. To this end, standard Newton-Raphson procedures have been employed for the solution of the resulting non-linear system of equations. The load-displacement curve for this particular node is shown in Figure 14.

In this case, the observed agreement between the PGD-ANM solution and that of the FE model is even higher than in the previous example. For the fourth order expansion, the agreement between the predicted load-displacement curves is almost exact.

Remember that our approach is based upon an off-line/online procedure such that, once the offline computation has been carried out, its solution is stored in the form of a series of 1D vectors that are evaluated in real time very efficiently. The solution thus computed is a multidimensional one, that is particularized online. This procedure is sketched in Figure 15.

CONCLUSIONS

The problem of PGD approaches to non-linear problems has been addressed here, with a particular emphasis on hyperelasticity. The development of suitable linearizations for complex non-linear problems formulated under a PGD framework has been a major focus of attention for our community in the last years. First, because consistent linearizations of the complex weak forms appearing in these formulations are far from being readily available. Previous works include explicit approaches to these problems [START_REF] Niroomandi | Real time simulation of biological soft tissues: a PGD approach[END_REF] but are difficult to generalize to arbitrary hyperelastic models or constitutive laws. The presented technique could be applied to any hyperelastic constitutive law, and, in view of existing results in the ANM, possibly to other non-linearities such as plasticity.

Here, we have proposed a combined PGD-ANM approach. The main advantages of this approach is that it produces a series approach to the solution that involves the same tangent operator for all the terms in the expansion. Thus, costly stiffness matrix updating procedures, which constitute nowadays to one of the main problems in the model order reduction community, are avoided. It has been shown how the proposed method provides with a very accurate solution to complex problems in the form of a metamodel that can be evaluated at very high feedback rates (even at real time rates on the order of 1 kHz) with very little computational cost. This is precisely one of the most attractive features of the just-proposed method. It combines an off-line/online procedure by which a truncated expansion of the multidimensional form of the displacement is obtained and stored in memory. After, these modes are evaluated in the online procedure of the method at extremely fast feedback rates. We believe that this philosophy opens new possibilities in the field of real-time simulation that deserve to be explored.

Other problems remain open, however. Notably, optimality of PGD approaches (i.e., under what circumstances a priori PGD modes are equivalent to a posteriori POD or SVD modes) is not well understood. It has been noticed in the examples throughout this work that, unlike those in previous explicit approaches, PGD-ANM modes seem not to be optimal. They are highly oscillating and numerous for a prescribed tolerance. The solution is simple, however. It suffices to obtain the SVD modes that show most of the energy of the system to obtain a very compact representation of the solution (the so-called projected PGD, see [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF]). But in any case, it seems pertinent to work in discerning what are the key ingredients of optimality in a PGD approach. This constitutes one of our current efforts of research.

In any case, although not optimal, the proposed method provides with a very competitive solution to highly demanding problems in applied sciences and engineering such as dynamic datadriven problems, real-time response even under very restrictive scenarios (haptic peripherals, for instance), simulation-driven control of structures and processes, and many others, where non-linear simulations are nowadays standard in industrial practice.

APPENDIX A: DERIVATION OF THE TANGENT STIFFNESS MATRIX FOR THE NEO-HOOKEAN CASE

For the neo-Hookean case explained in Section 5.2, the tangent stiffness matrix takes the following form:

K t D Z 0 B T DB C G T Q S 0 G Á d , ( A . 1 )
where

D D Â 1 J 2 0 C 1 0 M T 0 Ã C 2. ln J 0 / Â 1 J 2 0 C 1 0 M T 0 Á Q C 0 Ã (A.2)
now takes into account the material non-linearity and has a somewhat similar appearance to the Lagrangian elastic tensor at the initial state. J 0 and C 0 represent the Jacobian and right Cauchy-Green strain tensor of the initial solution. M 0 is obtained from the series expansion of the Jacobian, and contains minors of C 0 . Finally, Q C 0 is obtained from the series expansion of C 1 and contains components of C 0 , arranged in a particular way.

The geometrical non-linearities are included in the matrices B, G , and Q S 0 . B represents the usual strain-displacement matrix, G relates the nodal displacements u and the gradient of displacements vector, and, finally, Q S 0 represents a matrix that contains the initial stresses (we have chosen the same notation as in [START_REF] Cochelin | The asymptotic numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF]).

In the right-hand side of Equation ( 35), the non linear load vector f nl p is a vector containing information of material and geometrical non-linearities of all order problems ranging from order one to p 1. It can be written as Finally, it is necessary to expand Y D ln J and Z D J 2 by using Taylor series and the chain rule generalized to higher derivatives, At this point, a procedure entirely similar to that of Equation ( 36) is performed, leading to an equivalent expression of the displacement in terms of a power series in parameter a, where each term is composed by a finite sum of separable functions.

Y p D 1 2J 2

8 <:

 8 u mC1 .a/ S mC1 .a/ mC1 .a/

Figure 1 .

 1 Figure 1. Load-displacement curve (in terms of the load parameter ) for one particular load position of the pinched cube problem. Different solutions for different orders of expansion (p D 1, : : : , 5) compared with the FE solution by employing Newton-Raphson algorithms.

Figure 2 .

 2 Figure 2. Deformed configuration of the cube for one particular load position.

Figure 3 .Figure 4 .

 34 Figure 3. Deformed configuration of the beam for the load positioned at the beam tip. Note that no artificial gain in volume can be perceived.

Figure 5 .

 5 Figure 5. Geometry of the finite element model for the liver.

Figure 6 .

 6 Figure 6. Load-displacement curve (in terms of the load parameter ) for one particular load position of the liver palpation problem. Different solutions for different orders of expansion (p D 1, : : : , 6) compared with the FE solution by employing Newton-Raphson algorithms. Noticeably, the Newton-Raphson algorithm did not converge when solving the case D 1.

Figure 9 .

 9 Figure 9. Functions F 4 i .x/, for i D 1, 2, 3, and 20, respectively, for the fourth order expansion, in the simulation of the Kirchhoff-Saint Venant liver.

Figure 10 .

 10 Figure 10. Functions G 1 i .s/, for i D 1, 2, 3, and 20, respectively, for the first order expansion, in the simulation of the Kirchhoff-Saint Venant liver.

Figure 11 .

 11 Figure 11. Functions G 4 i .s/, for i D 1, 2, 3, and 20, respectively, for the fourth order expansion, in the simulation of the Kirchhoff-Saint Venant liver.

4 Figure 12 .Figure 13 .

 41213 Figure 12. Solution for the neo-Hookean beam bending problem at different expansion orders.

Figure 14 .

 14 Figure 14. Comparison between the PGD-ANM solution, at different expansion orders, and that for a standard finite element solution for a load at a particularized position. Neo-Hookean behavior.

Figure 15 .

 15 Figure 15. The result of the off-line PGD simulation is a multidimensional solution that depends on the load position. Thus, evaluating the solution for a prescribed load is entirely similar to performing a cut on a hyper volume N . In the figure, the red line on the blue liver represents approximately the path followed by the load over the surface of the liver. The displacement is magnified by a factor of 2.0.

3 ) 4 ) 1 p

 341 As in the stiffness matrix, S nlgeom p and S p represent the standard matrices found in literature when ANM is used to solve geometrical non-linear problems with linear materials. S nlmat p takes into account the material behavior,S nlmat p D . ln J 0 /In this equation, CC 0 represents the cofactor matrix of C 0 , and RCC p is a vector containing values of C ij of all problems from order one to p 1, obtained when the cofactor matrix of C is expanded in Taylor series,CC p D Q C 0 C p C RCC p ,RJ p is a summation of products of different components of C p and is obtained when the squared Jacobian is expanded in Taylor series, .J 2 / p D M T 0 C p C RJ p . ( A . 7 ) RS p collects terms concerning the expansion of Y D ln J and C 1 , collects terms concerning Z D J 2 and cofactor matrix of C expansions,

  Figure 7. Vertical displacement field of the Kirchhoff-Saint Venant liver for one particular position of the load.

		Y
		Z	X
		Uy
		-0.1
		-0.2
		-0.3
		-0.4
		-0.5
		-0.6
		-0.7
		-0.8
		-0.9
		-1
		-1.1
		-1.2
	Y		Y
	Z	X	Z	X
	Uy		Uy
	0.015		0.015
	0.01		0.01
	0.005		0.005
	0		0
	-0.005		-0.005
	-0.01		-0.01
	-0.015		-0.015
	-0.02		-0.02
	-0.025		-0.025
	-0.03		-0.03
	-0.035		-0.035
	-0.04		-0.04
	-0.045		-0.045
	-0.05		-0.05
	-0.055		-0.055
	-0.06		-0.06
	-0.065		-0.065
	(a)	(b)
	Y		Y
	Z	X	Z	X
	Uy		Uy
	0.035		0.022
	0.03		0.02
	0.025		0.018
	0.02		0.016
	0.015		0.014
	0.01		0.012
	0.005		0.01
	0		0.008
	-0.005		0.006
	-0.01		0.004
	-0.015		0.002
	-0.02		0
	-0.025		-0.002
	-0.03		-0.004
			-0.006
			-0.008
			-0.01
			-0.012
	(c)	(d)

Figure 8. Functions F 1 i .x/, for i D 1, 2, 3, and 20, respectively, for the first order expansion, in the simulation of the Kirchhoff-Saint Venant liver.

0

  .J 2 / p C RY p , andZ p D / 1 .J 2 / 2 , / 1 .J 2 / 2 ,

						1 0 J 4	.J 2 / p C RZ p ,	(A.10)
	where
	RY 1 D 0,
	RY 2 D	1 0 4J 4	.J 2 / 2 1 ,
	RY 3 D .J 2 RZ 1 D 0, 1 6J 6 0 .J 2 / 3 1 C 2 1 4J 4 0
	RZ 2 D	1 0 J 6	.J 2 / 2 1 ,
	RZ 3 D	1 J 8 0	.J 2 / 3 1 C 2	1 0 J 6	.J 2

: : :

ACKNOWLEDGEMENT

The authors would wish to thank the Spanish Ministry of Economy and Innovation (grant number, CICYT-DPI2011-27778-C02-01) for their support.