Sylvain Girbal
email: sylvain.girbal@thalesgroup.com

Jimmy Le Rhun
email: jimmy.lerhun@thalesgroup.com

BB-RTE: A BUDGET-BASED RUNTIME ENGINE FOR MIXED & TIME CRITICAL SYSTEMS

The safety critical industry is considering a shift from single-core COTS to multi-core COTS processors for safety and time critical computers in order to maximize performance while reducing costs.

In a domain where time predictability is a major concern due to the regulation standards, multi-core processors are introducing new sources of time variations due to the electronic competition when the software is accessing shared hardware resources, and characterized by timing interference.

The solutions proposed in the literature to deal with timing interference are all proposing a trade-off between performance efficiency, time predictability and intrusiveness in the software. Especially, none of them is able to fully exploit the multicore efficiency while allowing untouched, already-certified legacy software to run.

In this paper, we introduce and evaluate BB-RTE, a Budget-Based RunTime Engine for Mixed & Safety Critical Systems, that especially focuses on mixed critical systems. BB-RTE aims at guaranteeing the deadline of high-critical tasks 1) by computing for each shared hardware resource a budget in terms of extra accesses that the critical tasks can support before their runtime is significantly impacted; 2) by temporarily suspending low-critical tasks at runtime once this budget has been consumed.

I. INTRODUCTION

Safety-critical applications are usually characterized by stringent real-time constraints, making time predictability a major concern with regards to the regulation standards [START_REF]Functional safety of electrical, electronic, or programmable electronic safetyrelated systems[END_REF], [12], [START_REF]297: Software, electronic, integrated modular avionics (IMA) development guidance and certification considerations[END_REF] of the safety-critical and time-critical industries.

The industry is nowadays considering a shift from singlecore COTS (component off-the-shelf) to multi-core COTS processors for safety-critical and time-critical products. Such a shift is appealing both in terms of performance as well as in terms of size, weight and power (SWaP) [START_REF] Thomas | Lessons learned integrating COTS into systems[END_REF]. It also fits with the exponential growth in terms of performance requirements in the embedded domain.

However multi-core processor architectures are introducing new sources of time variations, and the solution providers can no longer rely on resource over-provisioning to enforce time predictability. As a consequence, the industry is facing a trade-off between performance and predictability [START_REF] Kirner | Obstacles in worst-case execution time analysis[END_REF], [START_REF] Mezzetti | On the industrial fitness of wcet analysis[END_REF]. As depicted in Figure 1, multi-core processors are characterized by shared hardware resources such as some levels of caches, the interconnect, the main memory or I/O controllers.

At hardware level, concurrent accesses on theses resources are arbitrated, introducing jitter at application level defined as timing interference [START_REF] Girbal | A complete toolchain for an interference-free deployment of avionic applications on multicore systems[END_REF]. Such interference, caused by electronic competition on shared hardware resources, are breaking the timing isolation principles required by the regulation standards [START_REF]Functional safety of electrical, electronic, or programmable electronic safetyrelated systems[END_REF], [12], [START_REF]297: Software, electronic, integrated modular avionics (IMA) development guidance and certification considerations[END_REF] of time-critical software.

Several papers [START_REF] Bin | Studying co-running avionic real-time applications on multi-core cots architectures[END_REF], [START_REF] Nowotsch | Leveraging multi-core computing architectures in avionics[END_REF] have studied the impact of timing interference on the Worst-Case Execution Time (WCET). Their authors have shown that not trying to tackle the interference problem leads to a performance loss for worst-case that can go far beyond the expected performance gain of using a multicore processor.

In a survey on Deterministic Platform Solutions [START_REF] Girbal | Deterministic Platform Software for hard real-time systems using multi-core COTS[END_REF] the authors presented various solutions allowing some level of time-deterministic usage of non-deterministic hardware platforms. Two families of solutions are presented: control solutions that aim at eliminating all interference by restricting the usage of the hardware platform, and regulation solutions that are focusing on keeping the impact of timing interference below a harmful level.

By introducing over-provisioning or complex runtimes, control solutions such as Marthy [START_REF] Jean | Ensuring robust partitioning in multicore platforms for ima systems[END_REF] or deterministic adaptive scheduling [START_REF] Fisher | Certifying Applications in a Multi-Core Environment: The World's First Multi-Core Certification to SIL 4[END_REF] fail to efficiently exploit the multi-core performance. Other solutions relying on execution models restrict too much the usage domain. For instance, The AER model [START_REF] Girbal | A complete toolchain for an interference-free deployment of avionic applications on multicore systems[END_REF] is easily applicable to distributed memory systems but lacks strong guarantees with regular shared memory systems.

Regulation solutions, on one hand, offer a better performance efficiency but lack the strong guarantees required by domains such as avionics. On the other hand they are perfectly adapted to less critical environments such as mixedcritical systems where high-critical tasks are running conjointly with low-critical tasks. In such systems, regulation solutions offer the possibility to degrade low-critical tasks to guarantee the timing behavior of high-critical tasks. To provide time guarantees, regulation solutions usually rely on budgeting the time with early deadlines or shifting time windows [START_REF] Kritikakou | Distributed run-time WCET controller for concurrent critical tasks in mixed-critical systems[END_REF].

In this paper, we present the Budget-Based RunTime Engine: a regulation solution relying on budgeting the number of shared hardware accesses to guarantee time properties. This is somehow similar to Memguard [START_REF] Yun | Memguard: Memory bandwidth reservation system for efficient performance isolation in multi-core platforms[END_REF], [START_REF] Sha | Single Core Equivalent Virtual Machines for Hard Real-Time Computing on Multicore Processors[END_REF] that preallocate application bandwidth requirement in term of memory accesses, but BB-RTE is especially focusing on the number of accesses, is extended to every shared hardware resource in the processor, and is performed per timeslot instead of as a whole.

To be able to do so, first we perform at design time an automatic standalone characterization of the critical applications to figure out the available budgets; and second we rely on this budget to take scheduling decision about non-critical applications at runtime. The paper is organized as follows: The principles of BB-RTE are presented in Section II and the full characterization and monitoring process is presented in Section III. In Section IV, we present the hardware target we considered, as well as the monitoring facilities we relied on. The mixed-critical software prototype is presented in Section V with its associated timing requirements. Finally, Section VI proposes an evaluation of the full process associated with the Budget-Based Runtime Engine, followed by a conclusion section.

II. BUDGET-BASE RUNTIME ENGINE: PRINCIPLES

The principles of the approach consist in determining, per timeslot, a maximum budget to allocate to low-critical applications in terms of resource accesses. When this budget is spent, low-critical applications are suspended until the next timeslot, as depicted in first timeslot of Figure 2.

Fig. 2. BB-RTE Principles

The budget is computed per timeslot and per hardware resource to ensure that the critical applications will be matching their deadline in each particular timeslot.

During the second timeslot of Figure 2, the low-critical applications manage to terminate without spending the whole allocated budget, again guaranteeing that the critical applications will match their deadlines during this timeslot.

With such a process based on budgeting, the major challenge consists in determining the budgets that guarantee the time behaviour of high-critical applications.

III. PROCESS TO DETERMINE AND ENSURE BUDGETING

The process to determine resource budgeting, depicted in Figure 3, is performed with two major steps.

First, both the architecture and the high-critical applications are characterized in an offline characterization step with regards to every hardware resources 1) to determine the total available resource budget; 2) to quantify the resource requirements of the high-critical task; 3) to figure out the maximum level of extra accesses to the resource before hampering the high-critical tasks.

Second, in a runtime regulation step, from the characterization information is inferred the maximum number of resource accesses allowed for non critical tasks, these tasks being suspended by a Budget-Based RealTime Engine once they reach this maximum number of total access during a given time slot.

A. Hardware characterization & budgeting

As introduced in Figure 3, in the first sub-step of the offline characterization phase, we perform a characterization of the target hardware platform. This hardware characterization consists in defining a set of low-level (assembly code) Stressing Benchmarks. Each of these stressing benchmarks is responsible for stressing a particular hardware resource of the selected multi-core target, by multiplying the number of accesses to this particular resource.

By progressively stressing each resource while monitoring both the execution time and the effective number of access to this resource thanks to Performance Monitor Counters (PMC) [START_REF] Sprunt | The basics of performance-monitoring hardware[END_REF], we are able to determine the maximum available bandwidth in terms of access to this resource, and that corresponds to the total available budget for the resource.

By iterating over all the potentially shared hardware resource, we obtain a vector of such total budgets that fully characterize the hardware limitations of the selected platform.

B. Critical application characterization & budgeting

During the second characterization sub-step of Figure 3, we are characterizing the usage made by high-critical applications of the shared hardware resources. To do so, we run the highcritical applications concurrently with the stressing benchmarks described above, progressively increasing the stressing level, and only monitoring the effect in terms of runtime of the high-critical applications.

It allows us to extract two different kinds of information:

First the required per-resource budget needed by our highcritical applications; and second, the level of extra resource access supported by our high-critical applications before being significantly slowed down. This supported extra accesses is the access budget that can be safely used by the low-critical applications.

The process to determine the acceptable level of slowdown and the associated extra access budget is depicted in Figure 4. The y-axis represents, during the current timeslot, the maximum observed runtime of the monitored application. The x-axis represents the extra access load performed on the associated hardware resource by the stressing benchmarks. The leftmost point in the chart corresponds to the application running alone in isolation (with a null stressing benchmark activity). It therefore corresponds to the classical WCET of the application running in isolation. The rightmost point in the chart corresponds to a permanent maximum load from the stressing benchmark actually preventing the monitored application to access the required resource.

Selecting an acceptable level of slowdown and the associated extra access budget is performed by selecting a point on the chart. The projection of this point on the yaxis then corresponds to the acceptable level of slowdown and the projection on the x-axis directly provides the extra access budget available for non-critical applications with an acceptable impact on critical applications.

We used two different techniques to select this point: First by directly selecting a maximum acceptable slowdown, but doing so for every hardware resource led us to too much overprovisioning, failing to fully exploit the multi-core efficiency. Second, we defined a maximum slope for the curve. This second solution allowed us to vary the level of slowdown relatively to the shape of the curve, focusing on the hardware resources the high-critical application was the most sensitive to.

By repeating this procedure for each shared hardware resource, we obtain again a budget in the form of a vector of both the required amount of access by the critical application (the leftmost point of the curve) and the number of supported extra access (the selected point of the curve).

Characterization steps involve large-scale experimentation due to the limited monitoring resources of multi-core processors. These architectures usually propose from tens to hundreds of hardware events that can be monitored, but only allow to monitor a few PMC at a given time (6 for ARM-v8 architectures for instance). As a consequence, testing all the countable PMC events involves many days of experimentation, and this characterizations phases have to be performed off-line. Anyhow this characterization only needs to be performed once per critical application for a particular hardware target, and such a characterization could also provide useful informations for the qualification or certification documents.

C. Budget-Based RunTime Engine (BB-RTE)

Once the characterization phases are over and all budgets have been gathered, both high-critical and low-critical applications can be deployed on the hardware target together with the runtime engine. As shown in Figure 3, the online regulation phase happens during the execution and consists in two substeps.

First, using the same process as the characterization phases, the low-critical applications are monitored with PMC counters, this time not to compute a budget, but to monitor the load in terms of hardware resource accesses.

Second, the BB-RTE runtime engine compares this load with the maximum extra budgets to decide if, during the current timeslot, non-critical task may continue executing or if they need to be suspended until the next timeslot. Doing so makes sure the slowdown of critical tasks does bot hamper their ability to match their deadlines.

One of the specific challenge of the BB-RTE runtime engine for time-critical systems is that time intrusiveness of the associated monitoring features has to be kept minimal, not to bias the time characterization results. Also the BB-RTE itself should make a minimal usage of shared hardware resources to not impact the resource access budgets. The final intrusiveness footprint of both the monitoring features and the BB-RTE engine will be presented in the result section.

IV. EXPERIMENTAL SETUP

The BB-RTE engine has been developed on top of PikeOS [START_REF] Sysgo | PikeOS 4.2: RTOS with hypervisor-functionality[END_REF] which is both a hypervisor and a real-time operating system relying on partitioning.

We relied on the METrICS toolsuite [START_REF] Girbal | METrICS: a measurement environment for multi-core time critical systems[END_REF] to perform the characterization steps and implemented the runtime engine as a native PikeOS partition, altering the scheduling in real-time as the monitoring information is gathered.

A. Target Architecture

The results presented in the Section VI of the paper were evaluated on an ARM Juno board [START_REF] Limited | Versatile Express: 64 Bit Juno r2 ARM Development Platform[END_REF] embedding a big.LITTLE architecture composed of a cluster of 2 high-performance Cortex A72 cores and a cluster of 4 more predictable Cortex A53 cores. The block diagram of the architecture is presented in Figure 5. The highlighted parts in the figure correspond to the shared memory path from the cores towards the main memory that are prone to timing interference. Other initiators like the MALI GPU, the PCI express links and the DMA controllers can also create contention on these shared resource, but for this paper we will only focus on memory contention from the cores. Alongside this memory path, the shared hardware resource are: The shared L2 cache that is shared by all the cores of each cluster. The CCI-400 cache-coherent interconnect (a high-bandwidth crossbar interconnect combining routing and coherency functions), and two DDR3 controllers served through the AMBA-compliant DMC-400 memory controller. The sizes of the memory structures alongside this path are presented in Table I. This information will be useful when trying to design stressing benchmarks dedicated at stressing a particular hardware resource.

B. Monitoring facilities

In Section III, we stated that the budget characterization is performed by using Performance Monitor Counters (PMCs).

The ARM-v8 architecture includes a Performance Monitors Unit (PMU), a non-invasive resource primarily used for debugging that provides information about the internal operations in the core. It includes a 64-bit cycle counter and 6 performance monitor counters able to count the occurrence of around 60 different events. Among these events we selected a subset of 13 events that correspond to either private or shared hardware resource composing the memory path, listed in Table II. The description appearing in the table is the level of details provided by the documentation on each countable event. A trial and error experimental process was therefore necessary to understand the exact meaning of each of these counters. inst_retired corresponds to the number of executed instructions whereas cpu_cycles corresponds to the number of CPU cycles measured so far. Therefore the ratio of the two correspond to the regular instruction per cycle (IPC) measurement commonly used to determine an application performance.

The events L1D_cache, L1I_cache and L2_cache indicate the number of accesses to the L1 data cache, L1 instruction cache and L2 cache by load/store requests respectively. The refill counterparts indicate how many time a cache line was obtained from a higher memory structure (e.g. the number of time a cache line was provided to the L1 cache by the L2 cache) so it is a good measure of the cache misses that indicates that a higher memory level in the datapath is accessed.

The write-back counterparts indicate how many time a modified cache line was written back to higher level memory so that the local cache line could be freed to fit a new data. This write back traffic, not appearing directly in the source code could be a significant part of the memory traffic.

The events bus_request and mem_request are respectively counting the requests on the local bus connecting core and caches, and the number of requests leaving the core to be send to the CCI-400 interconnect. This local bus traffic does not only correspond to memory accesses in the code, but also to instruction fetches and to the coherency traffic.

Finally, prefetch indicates the number of time a cache line of the L2 cache is filled due to an automatic hardware prefetch of the cache line. This feature has proven to be very problematic for our stressing benchmarks later presented in this section as the prefetching was limiting the stress level, as shown in the evaluation section of the paper. Some of the other hardware resource appearing in Figure 5 also provide some debugging support, but the documentation if often lacking or only available under NDA, and restricted to third-party providers of debugging probes such as Lauterbach. We definitely want to monitor these SoC-level events in the future, but we considered them out-of-scope for this paper.

C. Target Environment

In order to measure precise execution times and to sample hardware performance counters, we used our Measurement Environment for Time Critical Software (METrICS). Its main concept is to sample performance counters (including the cycle-count register) with a very short timing overhead, before and after the code sequences to monitor.

The METrICS toolsuite is composed of several elements: a kernel driver used to configure the hardware performance counters (as this requires privileged instructions), a library providing measurement probes to the application, and a Collector performing various initialization and transmission of measurement results. We heavily modified this latter component to include the part corresponding to the Run Time Engine. In addition, the METrICS suite includes host-side scripts and tools for measurement campaign automation, post-processing of raw data, and visualization.

The latency of a METrICS probe has been evaluated to be less than 392ns at worst. This whole probe thus has a comparable latency to a system call, and is precise enough considering the millisecond deadlines we have in the mixedcritical prototype presented in Section V.

The initial behaviour of the collector component was to run out of the monitored application operational cycles, to minimize the impact it had on the application timings. More precisely, it was 1) performing initialization and performance counter selection prior to running the monitored applications, 2) being completely suspended during the application operational cycles, 3) being activated again at the end of the application to dump the collected data to the host.

As appearing on Figure 6, we transformed in the BB-RTE context the collector component into the runtime engine. It is now also running during the application operational cycle to monitor resource usage, and to suspend low-critical tasks when they have consumed all their budgets. In the experimental context, we kept the third dumping phase to also collect the applications runtime and performance counter data including the number of deadline misses observed during the applicative phase. The application deployment also appear in Figure 6. We dedicated the two A72 cores to running both the operating system and the runtime engine. Several tests have shown that it was more efficient to run the RTE on the same core as the operating system due to the large number of system calls required to alter the scheduling. We therefore reserved the more deterministic A53 cores to the applications. Further deployment details for this cores will be provided in the evaluation section.

D. Stressing benchmarks

The mixed-critical applications we will use for the evaluation are described in Section V, but the characterization steps described in Figure 3 also involve stressing benchmarks. These benchmarks are simple applications performing repeated accesses to a shared resource. This involves executing a number of load or store instructions with regards to a memory region or a memory-mapped peripheral. In order to focus the stress applied by these benchmarks on a particular resource, we developed them in assembly language, thus allowing the greater control on the low-level behaviour of the stressing benchmark.

The accesses are parametrized in four ways: the direction of the transfer, the amount of data transferred, the address offset between consecutive accesses, and the repetition rate. The first parameter is either read or write. The second parameter, representing the buffer size, spans from 1KB to 2MB in power of 2 increments. This allows us to experiment the behaviour of data fitting or not in L1 and L2 caches. The address offset, called stride, is either 4 bytes for continuous accesses, 32 bytes to stress A72 L1 cache lines, or 64 bytes to stress A53 L1 and L2 cache lines. To vary the amount of accesses performed in a given time, we add a configurable number of NOP instructions (modeling calculations performed in isolation) per load or store instruction. This parameter spans from 0 (full stress) to 100 NOPs per access.

V. MIXED CRITICAL PROTOTYPE

To evaluate the Budget-Based runtime Engine, we set up a multi-domain mixed-critical software prototype composed of a high-critical application representative of the avionic domain, and low-critical control-command applications from the industry domain. The characterization of these applications as well as the evaluation of these applications running concurrently while supervised by the BB-RTE run-time engine will be provided in Section VI.

A. High-critical Application: Flyance

As the high-critical application, we selected Flyance, a mark-up Flight Management System (FMS) application mimicking the real-time behaviour of a regular FMS from the avionics domain.

The purpose of a Flight Management System (FMS) in modern avionics is to provide the crew with centralized control for the aircraft navigation sensors, computer-based flight planning, fuel management, radio navigation management, and geographical situation information. Taking charge of a wide variety of in-flight tasks, the FMS allows to reduce the workload of the flight crew.

The FMS is especially responsible for services that allow in-flight guidance of the plane. From pre-set flightplans (takeoff airport to landing airport), the FMS is responsible for plane localization, trajectory computation allowing the plane to follow the flightplan, and reaction to pilot directives. Flyance is composed of 25 time-critical tasks that are regrouped into different task groups as presented in Figure 7. The Sensors task group is in charge of generating all the localization data from various sensors. The Localization task group is in charge of computing the most probable position of the aircraft (BCP) by merging sensor information with different trustworthiness levels. The Nearest Airports task group continually builds a list of the nearest airports during the flight. The Flightplan task group is in charge of managing and processing modification requests on the flightplans that are preset routes used to guide the airplane. The Trajectory task group aims at computing both lateral and vertical profiles for the three flightplans set by the flightplan task. The lateral profile is composed of waypoints as well as leg information (path before, after and between the waypoints). The vertical profile provides altitude information (cruise altitude interceptions, crossing altitudes and slope angles) as well as performance information (estimated time of arrival, estimated fuel on board).

The FMS application also embeds a large Navigation Database that does not fit in any cache structure. It is both linearly and regularly accessed by a task from the Nearest airport task group, as well as randomly and sporadically accessed by tasks of the Flightplan task group. Accesses to this database in the main memory is very timing interference prone.

All the tasks composing the FMS have stringent real-time requirements presented in Table III. Additionally, aperiodic tasks have 100ms deadlines and can not be activated more than twice every 200ms. During the evaluation, one A53 core of the Juno board will be dedicated to running the Flyance application.

B. Low-critical Application: BiQuad

As the low-critical real-time application for the mixedcritical prototype, we selected Biquad: a control-command application implementing bi-quadratic on some streamed data acquired from analog sensors. Such an application is representative of classical distributed control system applications commonly found in the industry. The software architecture of the application is presented in Figure 8. BiQuad is composed of eight tasks: the Generator producing input data, either through generation on the first pass, or iterating on the received data on the next passes; the Splitter routes the data to the filtering tasks; the four Filters respectively apply their bi-quadratic filter on the data; and the Aggregator fuses the filtered data and sends it back as feedback to the generator task. The fused data is also periodically driving the Actuator task.

The timing requirement of this application is uniform across all the tasks: each task has a period of 200ms and a deadline equal to the next period.

For the evaluation, several instances of the BiQuad application can run at the same time, each running on a different A53 core.

VI. EVALUATION

This section first reports the budget characterization results obtained following the process described in Section III and summarized by Figure 3. it presents an evaluation of the runtime engine, by comparing the results in terms of deadline misses compared to an unregulated run of the same set of applications.

A. Platform characterization and Total available budget

Performing the hardware platform characterization of Figure 3 involves concurrently running stressing benchmarks until some saturation phenomenons are observed.

On the first A53 core, we run a monitored stressing benchmark with a fixed number of iterations and a full stressing profile (zero NOPs). On the other 3 A53 cores, we run the stressing benchmark as an infinite loop, starting it before the monitored one to be sure that the monitored benchmark is run under stress. For these stressing cores, we are varying the buffer size and the stride concurrently with the monitored core. We are also varying independently the number of NOPs (and therefore the amount of concurrent access to the resources).

By varying the stride, we impact which memory structure and therefore which part of the architecture will be exercised: With a 4-byte stride, we will will exercise the L1 cache maximizing the number of L1 hits, while strides larger than the L2 cache line size will mostly produce L1 and L2 cache misses, actually exercising the DDR controller and the memory.

Table IV shows the runtime results of such a characterization, varying the buffer size (how well the data will fit in the different cache structures), as the stressing level varying to maximum stress (no NOP instruction) to no stress (no load instructions). Selecting a 4-byte stride puts the L1 data cache under pressure, but this cache structure is private to each core. As a consequence, results presented in Table IV show very small variations while increasing the number of NOPs, the runtime only varying as the buffer size grows. Due to the high L1 hit ratio, having large buffer size not fitting in the L1 does not impacts the performance either.

Selecting larger stride values allows us to characterize further hardware components shared by all the A53 cores starting with the L2 cache and even with the A72 cores beyond the CCI 400.

As we were limited to count core-related events with the Performance Monitor Counters, we mainly focused on the level 2 cache, that is also the first hardware memory resource shared by all the cores of the A53 cluster. Whereas the results should allow us to observe timing interference with an expected runtime variability of fullystressed versus standalone larger than x4, we still observe very little difference between the two standalone and stressed versions.

Further tests, including unrolling the stressing benchmark loop to further increase the load or store ratio did not significantly change the above results. We therefore focused on the results corresponding to the maximum stress condition (2MB buffer size, with 0 NOP per iteration on the unrolled version) and studied the associated hardware counters presented in Table VI. In the worst case, out of the 35K data accesses to the L2 cache, 32K accesses are going to the external memory interface, meaning that our configuration successfully maximizes the number of L2 cache misses in order to maximize the interference on the interconnect. However, we can also observe 32K occurrence of prefetch, meaning that the hardware prefetcher successfully manages to capture the access pattern of our stressing benchmark and was able to anticipate nearly all of the external memory accesses. As a consequence, our stressing benchmark fails to effectively continuously stress the hardware resource.

In such a configuration the hardware prefetcher, instead of worsening runtime variability, smooths the contention on memory accesses. This seems to be a case of unexpected positive contribution to determinism from a dynamic, nondeterministic mechanism.

We tried to disable the hardware prefetcher, unfortunately we observed that it is periodically turned back on. We suspect the System Control Processor to be responsible for that. However, from our real time system, we have no control nor access on this particular core.

An alternative approach would be to develop another stressing benchmark, performing accesses that are more difficult for the prefetcher to predict. This would for example involve pseudo-random address increments.

B. Critical application characterization and Extra supported budget

The characterization of the critical application was performed by running the FMS application on the first A53 core. Again, thanks to METrICS, we instrumented the FMS application by inserting probes around each task composing the application, and at the level of a full operational cycle from the sensor task to the trajectory computation.

We first ran the the FMS application standalone, with the other core being idle, to collect the budget in terms of L2 accesses required by each timeslot of the application. Contrary to the stressing benchmarks, that are regular in their behavior and throughput requirements, the observed runtimes and resource requirements vary a lot with the FMS application, as all the tasks are not triggered in every 200ms timeslots due to heterogeneous periods.

The LOC C 4 task for instance is only executed in 1 out of 5 periods, whereas the LOC C 1 task is executed every period. As expected the timeslots corresponding to the least common multiple of all the periods, during which all the tasks must execute, is the one with the larger observed runtime and L2 access count, as pointed out by Figure 9. The next step, as defined in Figure 3, would have been to run the FMS application on the first A53 core, while running stressing benchmarks on the other A53 cores to build a figure looking like Figure 4. Such experiments led to the same issue as presented in the previous section, with the benchmarks failing to successively stress the architecture.

As a consequence, we failed to build such a figure with well identified asymtotes, and had to rely on an alternative technique:

Let R = [R 0 , R 1 , . . . , R N ] be the set of all the required access budget per timeslot for the critical FMS application running standalone without any activation of the asynchronous tasks. We set E = [E 0 , E 1 , . . . , E N ] the set of per-timeslot extra supported budget such as

E i = (max(R) -R i ) × 20%.
Doing so is not representative of the total budget available on the target platform, but this way the supported extra budget is inversely proportional to the resource usage performed by the critical application, as it should really be.

C. RTE Evaluation

As a baseline for the evaluation of the runtime engine we started to deploy the FMS application on the first A53 core, and the BiQuad applications on the remaining A53 cores. Doing so, we observed a deadline miss ratio (number of timeslot with a deadline miss compared to total number of timeslots) of 24.7%.

We then run the same application deployment, this time supervised by the Budget Based Runtime Engine. As a budget, we used the previously computed budget E. And again we ran the FMS application without any asynchronous tasks. The deadline miss ratio decreased to 2.6% while the number of slots with suspended non-critical tasks increased to 30.9%, proving the ability of the Runtime Engine to suspend lowcritical tasks when they are endangering the high-critical ones.

The high-critical application normally also encompass some asynchronous tasks that we ignored so far. However some of these sporadic task have a huge memory footprint, especially the tasks from the Trajectory task group that are performing random accesses to the 100MB large navigation database.

It would make no sense to have a preset static scenario, forcing in which timeslot each asynchronous task would be triggered. It would be similar to considering these tasks as periodic with a very large period. On the other hand a random scenario will raise a reproductability issue as we are gathering the timeslot statistics.

As a consequence, we created randomly two static scenarios for asynchronous task triggering, one was used during the characterization phases to compute a R extra supported budget set, and the second one was used while running under the supervision of the runtime engine. In such a scenario we observed a final miss ratio of 13.0% and a suspend ratio of 52.7%.

One one hand, the increase in suspend ratio is due to the over-provisioning of the timeslots with some asynchronous tasks in E , leading to unnecessary suspends at runtime. One the other hand, the increased deadline miss ratio could be explained by the lack of such provisioning in the timeslots were the asynchronous tasks finally occur.

D. Limitation of the approach

Beyond the issues related to the characterization phases because of the hardware prefetcher preventing our dedicated benchmarks to sufficently stress the hardware resources, we identified a set of issues and limitations while evaluating the runtime engine.

First relying on Performance Monitor Counters restricts us to the hardware resources within the cores. As a consequence we were not able to evaluate shared hardware resources such as the DDR controller. With the prefetchers disabled, furthermore increasing the stride of the stressing benchmarks would have caused memory page reloads at the level of the memory controller, causing extra delays. These sources of interference have not yet been evaluated.

A more fundamental limitation is the issue caused by the aperiodic tasks. When pre-computing budgets with characterization phases, we face the same issue as with machine learning with both over-learning, and outlier issues, and as a result are facing either deadline misses or unnecessary suspends. Also having per-timeslot budget is not practical for real applications running for hours. During our evaluation, the FMS ran for up to 5 minutes. With 200ms timeslots, that involves 1500 different budgets for this critical application. Also, because of unpredictable aperiodic tasks, it is not really possible to identify a shorter repeating patterns.

We performed test with an average timeslot budget, but doing so only reduces the number of deadline misses by 2.04%. Using a maximum timeslot budget on the other hand eliminates all the deadline misses but the low-critical tasks are systematically suspended.

VII. CONCLUSION AND FUTURE WORKS

In this paper, we presented a regulation solution based on budgeting that aims at guaranteeing the temporal behavior of high-critical tasks in a mixed critical system, while degrading the behaviour of non critical tasks.

The budgeting approach is performed offline and per timeslot. If the approach worked well with purely periodic tasks, high-overhead aperiodic or sporadic tasks caused some major problems: either over-provisioning if considered during the offline characterization phase, causing very low budget for high-critical tasks and as a consequence a high rate of lowcritical suspend rate; or causing high-critical deadline misses if occurring at runtime during an unprepared timeslot.

In single-core systems, it is common practice to have a dedicated periodic slot to deal with sporadic asynchronous tasks. The applicability of such a practice for multi-core architecture depends on the way the applications are deployed in the system to benefit from parallelism: an option is to parallelize inside applications / partitions, granting all the cores to a single application during each time slot. Such a deployment is compatible with the usual way of dealing with aperiodic tasks. But it also forces to re-write the applications, which performance will then be constrained by the Amdahl's law [START_REF] Amdahl | Validity of the single processor approach to achieving large scale computing capabilities[END_REF].

Another option is to run different independent applications, running partitions in parallel. If such a scheme is good for software development that could carry on producing singlethread applications, and even though it could bring better performance, exploiting the Gustafson's law [START_REF] Gustafson | Reevaluating amdahl's law[END_REF] and not being hampered by data dependency. It does not allow anymore to have dedicated timeslots for a particular kind of traffic or tasks, unless introducing costly synchronization.

As a consequence, before selecting the most adequate control or regulation solution to deal with timing interference on multi-core, a first step should be to consider the possible deployment of the applications, figuring out which kind of parallelism will be exploited. This choice has consequences on the timing interference level, on the visibility of interference (from white-box in case of intra-partition parallelism as they are coming from well known other tasks of the application, to black-box when coming from a potentially unknown independent application in case of inter-partition parallelism).

ACKNOWLEDGMENT

The research leading to this work has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 644080 (SA-FURE).

Fig. 1 .

 1 Fig. 1. Timing Interference in multi-core architectures

Fig. 4 .

 4 Fig. 4. Determining an acceptable level of slowdown and the associated extra access budget

Fig. 5 .

 5 Fig. 5. Simplified block diagram of the ARM Juno Board

Fig. 6 .

 6 Fig. 6. METrICS infrastructure in the context of BB-RTE

Fig. 8 .

 8 Fig. 8. Software architecture of the BiQuad application

Fig. 9 .

 9 Fig. 9. Variation of the runtime (in blue x) and L2 accesses (in red +) while running the FMS application standalone

TABLE I

 I

	.	JUNO BOARD MEMORY STRUCTURE SIZES
			A72 cluster A53 cluster
	Number of cores	2	4
	IL1 cache size	48KB	32KB
	IL1 cache line size	64B	64B
	IL1 associativity	3	2
	DL1 cache size	32KB	32KB
	DL1 cache line size	32B	64B
	DL1 associativity	2	4
	L2 cache size	2MB	1MB
	L2 cache line size	64	64
	L2 associativity	16	16

TABLE II

 II

	.	HARDWARE EVENTS MEASURED WITH PERFORMANCE
			MONITOR COUNTERS
	Performance counter Counting
	inst retired		Instruction architecturally executed
	cpu cycles		Cycles
	L1I cache		Level 1 instruction cache access
	L1I cache refill	Level 1 instruction cache refill
	L1D cache		Level 1 data cache access
	L1D cache refill	Level 1 data cache refill
	L1D cache wb	Level 1 data cache write-backs
	L2D cache		Level 2 cache access for data
	L2D cache refill	Level 2 cache refill for data
	L2D cache wb	Level 2 cache write-backs for data
	bus access		Bus access
	mem request	Memory access
	prefetch		Linefill because of prefetch

TABLE III

 III

	.	FLYANCE TIMING REQUIREMENTS
	Task	Period Deadline
	SENS C1	200 ms	200 ms
	LOC C1	200 ms	200 ms
	LOC C2	1600 ms 1600 ms
	LOC C3	1200 ms 1200 ms
	LOC C4	1000 ms 1000 ms
	TRAJ R1	200 ms	200 ms
	TRAJ R2	300 ms	200 ms
	TRAJ R3	300 ms	200 ms
	NEAR P 1 1000 ms 1000 ms

TABLE IV

 IV

	.	RUN-TIME VARIATION OBSERVED ON THE MONITORED
			CORE WITH A 4-BYTE STRIDE
					runtime (cycles)
	Buffer size	stress	min	25%	median	75%	max
	16384	none	12338	12352	12360	12369	13101
	16384	max	12338	12352	12360	12369	13175
	32768	none	24680	24786	24839	24959	26129
	32768	max	24679	24785	24838	24960	26272
	65536	none	49340	49951	50300	50338	52499
	65536	max	49339	49531	49752	50759	52795
	524288	none	393962	394096	394219	394633	417123
	524288	max	393978	394535	394785	395224	421488
	1048576	none	789994	791536	792055	792817	827872
	1048576	max	787995	791789	793538	795033	823978
	2097152	none	1599528	1601997	1603697	1606266	1751393
	2097152	max	1577078	1603203	1605720	1607560	1732193

TABLE V .

 V

		RUN-TIME VARIATION OBSERVED ON THE MONITORED
		CORE WHILE PERFORMING L2 CACHE MISSES
				runtime (cycles)
	Buffer size	stress	min	25%	median	75%	max
	16384	none	818	826	832	839	2883
	16384	max	818	824	831	839	4975
	32768	none	1728	1866	1902	1955	9739
	32768	max	1734	1883	1929	2006	7382
	65536	none	6772	8237	8885	10131	19378
	65536	max	6720	9080	9771	10777	20170
	524288	none	79625	107638	108660	111344	179419
	524288	max	97018	109665	110268	110897	177007
	1048576	none	169753	210845	214317	229595	366017
	1048576	max	207393	227506	237746	242661	353066
	2097152	none	478261	522073	526574	529720	725379
	2097152	max	500723	547663	551346	553397	705256
	Table V presents the running time distribution on the
	monitored core while varying again the buffer size and the
	number of NOP instructions, to compare standalone versus
	maximum stress deployments.

TABLE VI

 VI

	.	PERFORMANCE MONITOR COUNTERS RELATED TO THE
		L2 CACHE UNDER STRESSING CONDITION
		counter	min	median	max
		l2d cache	33645	34654	35129
		mem request 32473	32636	32863
		prefetch	30304	31716	32001


Anass Serhani

Ghislain Haine

Denis Matignon

Anisotropic heterogeneous n-D heat equation with boundary control and observation: II. Structure-preserving discretization

Keywords: Port-Hamiltonian Differential Algebraic System, Heat Equation, Structure Preserving Discretization, Partitionned Finite Element Method (PFEM), Boundary Control

 to the three associated port-Hamiltonian systems. Differential Algebraic Equations are obtained. The strategy proves very efficient to mimic the continuous Stokes-Dirac structure at the discrete level, and especially preserving the associated power balance. Anisotropic and heterogeneous 2D simulations are finally performed on the Lyapunov formulation to provide numerical evidence that this strategy proves very efficient for the accurate simulation of a boundary controlled and observed infinite-dimensional system.

INTRODUCTION

The port-Hamiltonian formalism provides a powerful way to model complex systems, including laws coming from different realms. It is based on the use of physically meaningful quantities, linked through two kind of equations, namely physical laws and constitutive relations. The system it then obtained by interconnecting each domain (mechanics, thermodynamics, electromagnetism, etc.) through their ports, describing the exchange of energy between them. Furthermore, the physical parameters are intrinsically taken into account in the model, allowing for anisotropy and heterogeneity. This formalism highlights an underlying mathematical structure, namely (Stokes-)Dirac structure, that guarantees a balance equation for the Hamiltonian functional: the power balance. The Hamiltonian often describes the physical energy, but other type of functionals can be used.

A recent topic of research is to provide accurate discretization methods (both in time and space) to preserve this powerful formalism, and especially the power balance. Several strategies have been proposed to this aim, as mixed finite element/Galerkin in [START_REF] Golo | Hamiltonian discretization of boundary control systems[END_REF]; [START_REF] Kotyczka | Structured discretization of the heat equation: Numerical properties and preservation of flatness[END_REF]; [START_REF] Kotyczka | Weak form of Stokes-Dirac structures and geometric discretization of port-Hamiltonian systems[END_REF], finite difference in [START_REF] Trenchant | Finite differences on staggered grids preserving the port-Hamiltonian structure with application to an acoustic duct[END_REF], finite volume in [START_REF] Kotyczka | Finite Volume Structure-Preserving Discretization of 1D Distributed-Parameter Port-Hamiltonian Systems[END_REF]; [START_REF] Serhani | Structure-Preserving Finite Volume Method for 2D Linear and Non-Linear Port-Hamiltonian Systems[END_REF]. For time-domain discretization, Celledoni and Høiseth (2017) closely examine the question.

In view of these works, it seems that the partitionned finite element method (PFEM), first introduced in Cardoso-Ribeiro et al. (2018), is the most accurate way to proceed in space, as this work will try to prove.

PFEM relies on the fact that Hamiltonian systems in infinite dimension are used to formalize distributed parameters systems, hence including partial differential equations (PDE) given by physical laws. The idea is to write a variational formulation of the pHs, and integrate by parts only one of the equations of the weak form, accordingly with the desired control (it has to "appear" in the boundary term of the integration by parts). Then we use finite element to obtain a finite dimensional pHs, with accurate properties with respect to the infinite-dimensional one. PFEM has been succesfully used on various type of physical systems, such as acoustic pressure waves or plates in Brugnoli et al. (2019b,c); Serhani et al. (2019d).

In the following, the heat equation is considered in a bounded open connected set Ω ⊂ R n , n ≥ 1, with mass density ρ(x), for all x ∈ Ω; -→ n denotes the outward unit normal. We suppose that this domain does not change over time, i.e. we work at constant volume. No chemical reactions are to be found in the domain. The quantities of interest are denoted as: u the internal energy density, -→ J Q the heat flux, T the local temperature, β := T -1 the reciprocal temperature, s the entropy density, -→ J S := β -→ J Q the entropy flux, and C V := du dT V the isochoric heat capacity.

In the companion paper Serhani et al. (2019a), the authors proposed to model the heat equation, controlled and observed at the boundary of the physical domain, by means of three different Hamiltonians. The first two are thermodynamically founded: the entropy and the internal energy, while the third one is less meaningful physically speaking, but enjoys an interesting behavior on a mathematical level. Reference to an equation (X) in this companion paper will be denoted by (I.X) in the present one.

We aim to apply the structure-preserving scheme PFEM to spatially discretize the three systems proposed in Serhani et al. (2019a) modeling the boundary controlled (and observed) heat equation. In the three cases, thermodynamical laws (i.e. constitutive relations) are needed to close the system. However, with the port-Hamiltonian formalism, the discretization of these relations can be postponed after the discretization of the physical laws. They need however to be carefully handled, suitably with the discretization of the port-Hamiltonian system.

The paper is organized as follows. In Section 2, PFEM is applied as a structure-preserving scheme on the three systems of Serhani et al. (2019a), and proves to be accurate for the corresponding Hamiltonians and its power balance. The Lyapunov case is then numerically experimented in Section 3 with 2D simulations to show the efficiency of the method, both on an isotropic and homogeneous case with analytically known solution and on an anisotropic and heterogeneous case. We end the paper with highlights and discussion about the pros and cons of PFEM compared to usual space discretizations.

A STRUCTURE-PRESERVING NUMERICAL METHOD: PFEM

This method is a recent structure-preserving method for pHs, first published in [START_REF] Cardoso-Ribeiro | A structurepreserving partitioned finite element method for the 2D wave equation[END_REF], and more recently developped in [START_REF] Cardoso-Ribeiro | A partitioned finite-element method (PFEM) for power-preserving discretization of open systems of conservation laws[END_REF].

The principle of the method lies on an integration by parts on a partition of the set of equations of the weak formulation of the continuous model, in order to make appear the control in the boundary integral.

Let us begin by a presentation of PFEM on the Lyapunov formulation. The method on the two other formulations will only be sketched, as most of the involved matrices will already be defined.

PFEM for the Lyapunov formulation

To begin with, let us write the weak formulation of (I.14

     Ω ρfuϕ dx = - Ω div -→ e Q ϕ dx, Ω -→ f Q • -→ ϕ dx = - Ω ---→ grad (eu) • -→ ϕ dx.
Integrating by parts the first line, one gets:

     Ω ρfuϕ dx = Ω -→ e Q • ---→ grad (ϕ) dx - ∂Ω -→ e Q • -→ n ϕ dγ, Ω -→ f Q • -→ ϕ dx = - Ω ---→ grad (eu) • -→ ϕ dx, (1)
leading to the choice of heat flux contol v

∂ := -- → e Q • - → n | ∂Ω .
Alternatively integrating by parts the second line

     Ω ρfuϕ dx = - Ω div -→ e Q ϕ dx, Ω -→ f Q • -→ ϕ dx = Ω eu div -→ ϕ dx - ∂Ω eu -→ ϕ • -→ n dγ, (2)
leads to the choice of temperature control v ∂ := e u| ∂Ω .

Assuming three finite-dimensional spaces

X := span{Φ} := span{(ϕ i ) 1≤i≤N }, X := span{ -→ Φ} := span{( -→ ϕ k ) 1≤k≤ -→ N }, X ∂ := span{Ψ} := span{(ψ m ∂ ) 1≤m≤N ∂ }, one can define f u,d (t, x) := Φ (x) • fu(t) := N i=1 f i u (t)ϕ i (x) fu(t, x), -→ f Q,d (t, x) := -→ Φ (x) • f Q (t) := -→ N k=1 f k Q (t) -→ ϕ k (x) -→ f Q (t, x),
and v ∂,d (t, γ) := Ψ (γ) • v ∂ (t) := N ∂ m=1 v m ∂ (t)ψ m (γ) v ∂ (t, γ).
Similarly, we define e u,d (t, x), -→ e Q,d (t, x) and y ∂,d (t, γ) in X , X , and X ∂ respectively.

The discrete weak formulation of (1) on

X × X × X ∂ reads    Mρfu(t) = De Q (t) + Bv ∂ (t), -→ M f Q (t) = -D eu(t), M ∂ y ∂ (t) = Ceu(t), (3)
with the following definitions

Mρ := Ω Φ • Φ ρ dx ∈ R N ×N , -→ M := Ω -→ Φ • -→ Φ dx ∈ R -→ N × -→ N , M ∂ := ∂Ω Ψ • Ψ dγ ∈ R N ∂ ×N ∂ , D := Ω ---→ grad (Φ) • -→ Φ dx ∈ R N × -→ N , B := ∂Ω Φ • Ψ dγ ∈ R N ×N ∂ , C := ∂Ω Ψ • Φ dγ ∈ R N ∂ ×N ,
where ---→ grad (Φ) means the rectangular matrix of size N ×n, constituted of the gradient of each functions of Φ, i.e.

---→ grad (Φ) := ---→ grad ϕ 1 , • • • , ---→ grad ϕ N .
Note that the mass matrices M ρ , -→ M and M ∂ are obviously symmetric. Furthermore, M ρ discretizes the inner product of L 2 ρ (Ω). Finally, C = B , hence (3) can be rewritten under the classical finite-dimensional pHs form, including the interconnection ports into the structure

M d -→ f d = J d -→ e d , (4)
with M d a symmetric matrix and J d a skew-symmetric matrix, both of size

N + -→ N + N ∂ . More precisely -→ f d := fu, f Q , -y ∂ , -→ e d := eu, e Q , v ∂ , M d := Mρ 0 0 0 -→ M 0 0 0 M ∂ , J d := 0 D B -D 0 0 -B 0 0 .
Note that the skew-symmetry of J d and (4) imply

-→ e d M d -→ f d = 0. ( 5
)
Now we want to define the discrete Hamiltonian and provide a discrete version of the balance equation (I.15).

Since

f u (t, x) := ∂ t u(t, x) in L 2 ρ (Ω), this leads to Ω ρ(x)fu(t, x)ϕ(x) dx = Ω ρ(x)∂tu(t, x)ϕ(x) dx.
With the discrete family Φ, this becomes

Mρfu(t) = Mρdtu(t).
(6)

In the same way, e u := u(t,x) C V (t,x) in L 2 ρ (Ω) leads to

Mρeu(t) = M ρC -1 V (t)u(t), (7)
where M ρC -1

V (t) := Ω Φ • Φ ρ C V (t,•) dx.
Denoting u d (t, x) the discrete counterpart of the internal energy density in X , one can define the discrete Hamiltonian by

H d (t) := 1 2 Ω ρ(x) (u d (t,x)) 2 C V (t,x) dx, leading to H d (t) = 1 2 u (t)M ρC -1 V (t)u(t).
Hence

dtH d (t) = dtu (t)M ρC -1 V (t)u(t) + 1 2 u (t)dtM ρC -1 V (t)u(t).
From ( 6) and ( 7), one immediatly gets

dtH d (t) = fu (t)Mρeu(t) + 1 2 u (t)dtM ρC -1 V (t)u(t).
Since dtM ρC -1

V (t) = -M ρC -1 V (t)dtM -1 ρC -1 V (t)M ρC -1 V (t)
, and ( 7)

gives u (t)dtM ρC -1 V (t)u(t) = -eu (t)MρdtM -1 ρC -1 V (t)Mρeu(t).
The two above equalities and (5) lead to the discrete counterpart of (I.15):

dtH d (t) = -e Q (t) -→ M f Q (t) + v ∂ (t)M ∂ y ∂ (t) - 1 2 eu (t)MρdtM -1 ρC -1 V (t)Mρeu(t). (8)
At this stage, constitutive relations have not been used. Assuming C V to be time-invariant and u = C V T , as in the Dulong-Petit model, leads in a weak sense to

Ω ρ(x)u(t, x)ϕ(x) dx = Ω ρ(x)C V (x)T (t, x)ϕ(x) dx, giving Mρfu(t) = M ρC V dtT (t), (9)
where

M ρC V := Ω Φ • Φ ρC V dx.
Of course, the definition of e u given by ( 7) then simply becomes

Mρeu(t) = MρT (t). ( 10
)
Fourier's law (I.5) writes

Ω -→ J Q (t, x) • -→ ϕ (x) dx = - Ω λ(x) • ---→ grad (T (t, x)) • -→ ϕ (x) dx, leading to -→ M e Q (t) = -→ Λ f Q (t), (11)
where

-→ Λ := Ω -→ Φ • λ • -→ Φ dx ∈ R -→ N × -→
N , which is symmetric by assumption on λ. Hence (8) becomes

dtH d (t) = -f Q (t) -→ Λ f Q (t) + v ∂ (t)M ∂ y ∂ (t), (12)
the discrete counterpart of (I.16).

Following the same strategy, we obtain the discrete formu-

lation of (2) in X × X × X ∂    Mρfu(t) = De Q (t), -→ M f Q (t) = -D eu(t) + Bv ∂ (t), M ∂ y ∂ (t) = Ce Q (t), (13)
with the following additional definitions

D := - Ω Φ • div -→ Φ dx ∈ R -→ N ×N , C := B := - ∂Ω -→ Φ • -→ n • Ψ dγ ∈ R -→ N ×N ∂ ,
where div

- → Φ (resp. - → Φ • - → n ) is the (column) vector of size -→
N , constituted of the divergence (resp. the normal trace) of each functions of the family -→ Φ.

Note that it is again possible to write this system under the form (4), with

J d :=   0 D 0 -D 0 B 0 -B 0   .

PFEM for the energy formulation

Now the aim is to discretize the system related to the internal energy U as Hamiltonian, i.e. system (I.11). Its variational formulation reads

             Ω ρfsϕ dx = - Ω div -→ e S ϕ dx - Ω eσϕ dx, Ω -→ f S • -→ ϕ dx = - Ω ---→ grad (es) • -→ ϕ dx, Ω fσϕ dx = Ω esϕ dx.
Practically, it seems legitimate here to control the temperature, and thus integrate by parts the second line

             Ω ρfsϕ dx = - Ω div -→ e S ϕ dx - Ω eσϕ dx, Ω -→ f S • -→ ϕ dx = Ω es • div -→ ϕ dx - ∂Ω es -→ ϕ • -→ n dγ, Ω fσϕ dx = Ω esϕ dx.
Following the same strategy as for the Lyapunov formulation, the following matrix form is computed

     Mρfs(t) = De S (t) -M eσ(t), -→ M f S (t) = -D es(t) + Bv ∂ (t), M fσ(t) = M es(t), M ∂ y ∂ (t) = Ce S (t),
where

M := Ω Φ • Φ dx ∈ R N ×N . ( 14
)
The pHs form (4) then reads

   Mρ 0 0 0 0 -→ M 0 0 0 0 M 0 0 0 0 M ∂       fs f S fσ -y ∂    =    0 D -M 0 -D 0 0 B M 0 0 0 0 -B 0 0       es e S eσ v ∂    . The discrete Hamiltonian is U d (t) := Ω ρ(x)u d (t, x) dx = u (t)M ρ 1, where 1 ∈ R N is a vector full of ones.
The aim is to show that U d enjoys a similar balance equation as U , i.e. similar to (I.13).

Gibbs formula (I.3) in a weak form reads

Ω ρ(x)∂tu(t, x) dx = Ω ρ(x)T (t, x)∂ts(t, x) dx,
which gives at a discrete level with the pHs variables

dtu (t)Mρ1 = fs (t)Mρes(t). ( 15
)
From the definition σ :=

---→ grad (β) • -→ J Q , one easily gets that T σ = - ---→ grad (T ) • -→ J S . This leads to Ω T (t, x)σ(t, x) dx = - Ω ---→ grad (T (t, x)) • -→ J S (t, x) dx,
which then reads once discretized, in the pHs variables -eσ (t)M fσ(t) = e S (t) -→ M f S (t). This implies that (5) reduces to es (t)Mρfs(t) = v ∂ (t)M ∂ y ∂ (t). From (15), one gets the discrete counterpart of (I.13), namely:

dtH d (t) = v ∂ (t)M ∂ y ∂ (t).
To close the system, and be able to solve it, we need accurate constitutive relations, according to the physical system under study, as done in the Lyapunov case.

PFEM for the entropy formulation

Now, we discretize the system related to the entropy S as Hamiltonian, i.e. system (I.8). Practically, it seems legitimate here to control the boundary heat flux. In pHs variables, this is exactly the same system to discretize as for the Lyapunov formulation with boundary heat flux control. The difference lies in the definition of e u and -→ f Q . Thus we have again system (13).

The discrete Hamiltonian is S d (t) := Ω ρ(x)s d (t, x) dx = s (t)Mρ1. From Gibbs formula (I.3), as for ( 15), we get dts (t)Mρ1 = fu (t)Mρeu(t). From ( 5):

eu (t)Mρfu(t) + e Q (t) -→ M f Q (t) = v ∂ (t)M ∂ y ∂ (t). Finally, the definition σ := ---→ grad (β) • -→ J Q can be directly discretized by σ (t)M 1 = -f Q (t) -→ M e Q (t).
Hence the discretized counterpart of (I.9)

dtS d (t) = v ∂ (t)M ∂ y ∂ (t) -f Q (t) -→ M e Q (t).

SIMULATION RESULTS

We now test PFEM on 2D examples. Simulations are performed on four cases: with heat flux boundary control or temperature boundary control, both for an isotropic homogeneous case where an analytical solution is known, and an anistropic heterogeneous case. We work with the Lyapunov Hamiltonian H under the hypotheses leading to (12). This choice was made as this is the usual "energy functional" in the applied mathematics literature.

The chosen geometry is the rectangle Ω := (0, 2) × (0, 1).

We propose to use P 2 Lagrange finite element for Φ, RT 1 Raviart-Thomas finite element for -→ Φ and P 2 Lagrange finite element for Ψ. The corresponding degrees of freedom (Dof) are represented on Figure 1. This leads to N = 204 Dof for Φ, -→ N = 470 Dof for Φ and N ∂ = 50 Dof for Ψ. The temperature T is approximated in the basis Φ. The discretization process will lead, either to an ordinary differential equation (ODE), or to a differential algebraic equation (DAE). To perform the simulations, we also need solvers to integrate in time these finite-dimensional differential equations. We use existing libraries to do it, without looking for time schemes which do not add numerical dissipation (as symplectic schemes for ODEs).

The simulations are carried out in Python 3, using FEniCS [START_REF] Alnaes | The FEniCS Project Version 1.5[END_REF]) for the finite element libraries. The time integration for the ODEs is performed with the Runge-Kutta RK45 scheme, and with IDA (SUNDIALS) for DAEs. Note that the resolution of the DAE needs compatible initial data, which can be a drawback.

Analytical solution

In this case, all the parameters are taken constant equal to one, i.e. ρ = C V ≡ 1, λ ≡ 1 0 0 1 . Let x := x 1 , x 2 and

T (t, x) := 4t + x 2 + 3x 1 -5x 2 . Moreover, define -→ J Q (t, x) := - ---→ grad (T (t, x)) = -2x 1 + 3, 2x 2 -5 . It is clear that ∂tT (t, x) = -div -→ J Q (t, x) = ∆T (t, x), ∀t ≥ 0, x ∈ R 2 .
Let us define the following function as boundary control:

v ∂ (t, γ) :=    5, x 2 = 0, 7, x 1 = 2, -3, x 2 = 1, -3, x 1 = 0, ∀t ≥ 0, γ ∈ ∂Ω. (16)
Then T defined above is the solution of the following heat equation with boundary heat flux control

∂tT (t, x) = ∆T (t, x), ∀t ≥ 0, x ∈ Ω, T (0, x) = x 2 + 3x 1 -5x 2 , ∀x ∈ Ω, ---→ grad (T (t, γ)) • -→ n (γ) = v ∂ (t, γ), ∀t ≥ 0, γ ∈ ∂Ω. (17)
The associated observation is then given by

y ∂ (t, γ) =      4t + x 2 1 + 3x 1 , x 2 = 0, 4t + x 2 2 -5x 2 + 10, x 1 = 2, 4t + x 2 1 + 3x 1 -4, x 2 = 1, 4t + x 2 2 -5x 2 , x 1 = 0, ∀t ≥ 0, γ ∈ ∂Ω. (18)
A straightforward computation gives

H(t) = 16t 2 + 52 3 t + 1301 90 . Hence dtH(t) = 32t + 52 3 . Also Ω ---→ grad (T (t, x)) 2 dx = 256
3 , and ∂Ω v ∂ (t, γ)y ∂ (t, γ) dγ = 32t + 320 3 , from which one verifies (I.16).

Boundary heat flux control

Equation ( 17) is the prototype of linear parabolic PDE with Neumann boundary condition. Discretizing it by the finite element method is a classical exercise which can be found in any lecture notes on FEM. It only makes use of the family Φ. It gives rise to the following ordinary differential equation (ODE)

M dtT (t) = -AT (t) + L(t), ∀t ≥ 0,
where M is defined in (14),

A := Ω ---→ grad (Φ) • ---→ grad (Φ) dx ∈ R N ×N ,
and

L(t) := ∂Ω v ∂ (t, γ)Φ (γ) dγ ∈ R N .
This gives a first way to simulate T , called ODE-FEM in the sequel.

The second way is to solve the DAE constituted of (3)-( 9)-( 10)-( 11). This method will be refer to as DAE-PFEM.

Finally, after substitutions, the above DAE leads also to the following ODE for the temperature

M ρC V dtT (t) = -D -→ M -1 -→ Λ -→ M -1 D T (t) + Bv ∂ (t),
i.e. using isotropy, homogeneity and the values given to the parameters: M dtT (t) = -D -→

M -1 D T (t) + Bv ∂ (t). We will refer to this method as ODE-PFEM.

Remark that in ODE-PFEM, the matrix D -→ M -1 D is not sparse, contrarily to matrix A in ODE-FEM. One can appreciate on Figure 2 the efficiency of the two PFEM approaches, leading to a relative error of 10 -4 , and still decreasing in time, whereas ODE-FEM seems to grow asymptotically to a relative error close to 10 -1 .

Furthermore, it has to be noted that ODE-FEM and ODE-PFEM have the same size and are integrated in time with the same scheme. Thus, despite the non-sparsity of the finite element matrix D -→ M -1 D leading to an overhead in CPU time and memory storage, there is still some interest in following this approach.

Boundary temperature control Now, the control v ∂ and the observation y ∂ respectively defined in ( 16) and ( 18) are switched, i.e. the heat equation is controlled by a prescribed temperature at the boundary of the domain (y ∂ ) and the heat flux is observed (u ∂ ).

The ODE-FEM approach has now a drawback, since the boundary temperature does not show up when integrating by part the weak formulation of (I.6). Thus, one needs a suitable scheme to impose the boundary control, such as lifting (time consuming since the control is timedependent), or Lagrange multiplier. We use the latter in our simulations.

The DAE-PFEM approach only consists of solving the DAE constituted of ( 13)-( 9)-( 10)-( 11), while substitutions allows again an ODE-PFEM approach

M dtT (t) = -D -→ M -1 D T (t) + D -→ M -1 Bv ∂ (t).
Figure 3 shows again the efficiency of PFEM, but only for the DAE resolution which gives a relative error at the computer precision. Remark that the error for ODE-PFEM, despite is value (10 -13 ) is not better than ODE-FEM. However, in long time, ODE-PFEM also reaches the computer precision while ODE-FEM stands still.

Anisotropic heterogeneous case

For this case, the parameters are taken to be

ρ(x) := x(2 -x) + 1, C V := 3, λ(x) := 5 + x 1 x 2 (x 1 -x 2 ) 2 (x 1 -x 2 ) 2 3 + x 2 x 1 + 1 .
The initial temperature is taken as a centered Gaussian, such that both temperature and heat flux are null at the boundary. Then a control is imposed either for the heat flux, or the temperature. It is taken to be under the form v ∂ (t, γ) = f (t)g(γ), with f (t) ∼ t t+1 and g(x 1 , x 2 ) = x 1 + x 2 at the boundary.

As an analytical solution would be difficult, not to say impossible, to compute, we compare the L 2 (∂Ω) relative error between the target control, and the effective value at the boundary of the approximation. Boundary temperature control Figure 5 shows that the control is far better imposed with the classical Lagrange multiplier method. This can be explained by the fact that the Lagrange multiplier method strongly imposes the condition, while PFEM approaches weakly impose it.

CONCLUSION AND PERSPECTIVES

To conclude this paper, we summarize what has been done. The heat equation has been written in the pHs formalism, with three different Hamiltonian (Serhani et al. (2019a)). Two of them are thermodynamically founded, namely the entropy and the internal energy, and the third is a weighted L 2 norm, usually used in applied mathematic. PFEM has been applied on the obtained systems and it proves to be an accurate way to spatially discretize them, leading to finite-dimensional pHs. In particular, this leads to an accurate discretization of the chosen Hamiltonian. The main difficulty lies in the discretization of the constitutive relations, which has to be suitably adapted to the discretization of the flows and efforts. We finally perform simulations in 2D, under usual hypothesis (Fourier's law and time-invariant isochoric heat capacity) for the Lyapunov formulation. We show on the one hand, the efficiency of PFEM to simulate the system, first compared to an analytical solution in a simple case, and for an anisotropic heterogeneous case. On the other hand, we observed that solving directly the DAE system seems to be more relevant, although it involves a bigger system, and suitable time schemes.

Some numerical analysis should be made, as this is done in [START_REF] Serhani | Numerical analysis of a semi-discretization scheme for the n-D wave equation in port-Hamiltonian system formalism[END_REF] for the wave equation, leading to an optimal choice for the finite element families.

Furthermore, a deeper investigation of the theoretical relations between the finite-dimensional pHs and the continuous one would be of great benefit, in particular concerning the family of Dirac structures (as function of the mesh size parameter), for a better understanding of the underlying geometrical structures; see e.g. Serhani et al. (2019c).

Figure 1 .

 1 Figure 1. The used mesh, and the relative Dof.

Figure 2 .

 2 Figure 2. Relative errors between the discrete Hamiltonians and the analytical one (heat flux control).

Figure 3 .

 3 Figure 3. Relative errors between the discrete Hamiltonians and the analytical one (temperature control).

Figure 4 .

 4 Figure 4. Relative errors in L 2 (∂Ω) between the target and effective boundary heat flux control.Boundary heat flux control It can be seen on Figure4, that DAE-PFEM has the best behavior. This can easily be explained. Since ODE-FEM and ODE-PFEM only approximates the temperature T , a post-processing is needed to compute the boundary heat flux, while it is directly accessible with DAE-PFEM.

Figure 5 .

 5 Figure 5. Relative errors in L 2 (∂Ω) between the target and effective boundary temperature control.

This work has been performed in the frame of the Collaborative Research DFG and ANR project INFIDHEM n • ANR-16-CE92-0028 (http://websites.isae.fr/infidhem).