A Mensch

F Charpillet

Scheduling in the REAKT kernel : combining predictable and unbounded computations for maximising solution quality in real-time knowledge-based systems

Keywords:

La complexité croissante des applications temps réel a renouvelé l'intérêt porté aux techniques de l'Intelligence Artificielle dans ce domaine. Cet article présente le Noyau d'Exécution REAKT, un environnement logiciel adapté aux besoins des applications temps réel à base de connaissances, et en particulier les techniques d'ordonnancement utilisées pour faciliter l'exécution conjointe de traitements déterministes et non déterministes tout en garantissant les temps de réponse.

Les techniques utilisées reposent sur le principe de raisonnement progressif, une approche qui permet à une application d'obtenir rapidement une première solution, puis de l'affiner tant que le temps disponible le permet. Un modèle de tâches est proposé pour permettre l'intégration aisée d'actions reflexes et de raisonnements complexes. Un algorithme permettant de maximiser le temps disponible pour les traitements complexes non déterministes tout en garantissant les temps de réponse des tâches temps réel. Cet algorithme fournit une réponse au problème d'ordonnancement joint de tâches périodiques, sporadiques et optionnelles.

THE NEED FOR MORE " INTELLIGENT " REAL-TIME SYSTEMS

Over the last few years, the usefulness and maturity of Artificial Intelligence technologies, and in particular of knowledge-based systems, have been demonstrated by the ever growing number of real-life applications using them. In the domain of real-time systems (process control, communication network management, satellites and avionics, patient monitoring, C 3 I systems, robotics, data analysis...), these techniques have been considered with a renewed interest, as they appear as a promising approach to cope with the increasing complexity of the systems to be controlled [Laffey et al. 88].

Although covering a wide range of domains, complex real-time applications show several common features :

• They are tightly coupled with a physical system (plant, network, satellite, robot, human body, battlefield...) which has to be controlled. This system evolves dynamically under the influence of external constraints which are not completely known by the control application. The behaviour of such an application is therefore strongly directed by the evolution of the controlled system, which generates data throughput and response times requirements. • They run in an environment which presents characteristics difficult to handle for classical real-time systems. This environment is often very dynamic, as in the case of communication networks or industrial processes, ill-modelled, as in the case of patient monitoring or financial analysis systems, or ill-defined, as in the case of autonomous vehicles or satellites. Such an environment increases the complexity of the control applications, as it requires from the application the ability to adapt itself to the evolving environment and to handle correctly time-varying, fuzzy, uncertain and incomplete information. • They have to provide results which go beyond the capabilities of today's real-time systems. Although current applications perform low-level tasks (such as data acquisition, data filtering and processing, detection of abnormal behaviour) very well, they usually depend on a human operator for higher-level tasks (such as diagnosis, situation assessment, action planning), as those tasks require a deep understanding of the mechanisms involved in the evolution of the controlled system. This approach becomes increasingly difficult with the growing complexity of the controlled systems, for practical, safety and economical reasons. These complex applications thus require the integration of advanced technologies in classical realtime systems in order to better assist, and sometimes even replace, the human operator in the decision process.

We present in this paper some results of the REAKT and REAKT II ESPRIT projects (EP 5146 and 7805). The primary objective of those projects was to develop a set of tools and the associated methodology for applying knowledge-based systems in realtime domains. The projects were to produce definitions, specifications and prototypes of various techniques, to be eventually integrated into a toolkit for the efficient development, deployment and maintenance of knowledge-based modules that can be embedded into real-time applications.

Section 2 presents an overview of the REAKT Real-Time Kernel, a software architecture well-adapted to the execution requirements of real-time knowledge-based systems. Section 3 describes the model of application activities which is used for scheduling tasks in the kernel. Sections 4 to 7 then detail the algorithms which support the joint scheduling of predictable real-time tasks and unbounded optional computations. Finally, Section 8 presents some conclusions on advantages of the approach.

THE REAKT KERNEL

REAKT is a software environment for developing and delivering real-time knowledgebased systems (RTKBSs). The REAKT environment provides both tools and a complete methodology to ensure the efficient development of applications with the requirements presented in the previous section. This environment is made of three main components:

• REAKT gets its unique real-time capabilities from its Real-Time Kernel, which provides all the necessary techniques for co-operative and real-time problem solving. This paper focuses on the REAKT Real-Time Kernel. More information on the REAKT Development Toolbox and the REAKT Application Methodology can be found respectively in [Galan et al. 93] and [Fjellheim et al. 94].

The REAKT Kernel [Mensch et al. 94] The Blackboard approach is widely recognised as a powerful framework which allows:

• the scaling limitations of expert systems to be overcome through a modular approach; • heterogeneous knowledge representation formalisms to be integrated;

• the application development life cycle and runtime execution to be better controlled.

While providing the benefits of the Blackboard approach, REAKT extends the classical model along several directions to provide an architecture adapted to real-time operation.

A REAKT application is composed of a set of tasks (Figure 2.1), each task being an independent thread of control in charge of executing one or several user-defined knowledge sources to contribute to the global solution of a problem. One or several input tasks can also be defined to perform data acquisition. Tasks can be either permanent or created dynamically in response to events. Tasks can be associated with deadlines, which represent the time at which the solution provided by a task must be available to avoid a system failure.

Figure 2.1 -Architecture of the REAKT Real-Time Kernel

REAKT uses a control strategy for task scheduling based on deliberative scheduling techniques [Garvey and Lesser 94] that offer efficient mechanisms to trade off computation time for quality of result. The proposed model follows the progressive reasoning paradigm [Charpillet and Boyer 94]. Progressive reasoning allows a system to obtain a first solution quickly, and then refine it while time is available. Two main types of tasks can be defined in a REAKT application:

• Real-time tasks, which must be permanent tasks with a bounded worst-case execution time and a task activation pattern which is either periodic or features minimal interarrival times. With these hypotheses, the REAKT scheduler is able to determine off-line whether a schedule meeting all real-time task deadlines is feasible, thus guaranteeing all hard deadlines in an application. Typical activities which can be performed at the first level include data acquisition, recognition of abnormal situations, reflex actions, output,... • Optional tasks (called intentions in Figure 2.1), which improve the solution quality while time is available, using the opportunistic control mechanism of Blackboard systems. Optional tasks are either refinements of real-time activities, or dynamically created tasks with only soft deadlines. These tasks are scheduled according to their deadlines and the quality of the solution already obtained, but their execution cannot be guaranteed, as neither their number nor their execution time need to be bounded.

THE REAKT TASK MODEL

REAKT provides a task model to describe the activities of an application. The introduction of a task model is required for two main reasons:

• No guarantee on the response times can be provided if the worst-case event arrival pattern in the application cannot be precisely estimated, i.e. the maximum load of the system is unknown. • The progressive reasoning approach requires related real-time and optional tasks to be declared to the system.

In the REAKT task model, each task τ is decomposed into three parts: mandatory, optional and action (Figure 3.1) [Audsley et al. 91]. The mandatory part, activated at Rm τ , is in charge of providing a first-level solution with a guaranteed response time; the optional part then tries to improve this result, using more complex reasoning mechanisms; the action part is activated before the task deadline D τ , possibly preempting the optional part, and uses the best available solution to execute the necessary actions.

Figure 3.1 -Decomposition of a task into mandatory, optional and action parts

Important characteristics of a task include:

• The type of the task: either periodic, when the task activation pattern repeats itself, or sporadic, when the task is triggered once in response to a particular event.

• The worst-case execution time of both the mandatory and the action parts of the task.

• The temporal characteristics of the task: period and deadline for a periodic task, minimum interarrival time and deadline for a sporadic task. • The name of the intention structure to be used to execute the optional part. This structure is described in [Lalanda et al. 92].

The REAKT Kernel relies on a two-layer control strategy, with a first level in charge of scheduling the real-time mandatory and action parts of each task, and a second level responsible for optimising the solution quality by running optional parts in the remaining time. The remaining sections of this paper focus on the algorithms used in the real-time control layer.

THE REAL-TIME CONTROL LAYER

The role of the real-time control layer is twofold: it must guarantee the response time of the mandatory and action parts of all tasks in an application, but should also maximise the amount of time available for optional second-level activities, in charge of improving the solution quality. For a given task, the basic principle is to first execute the mandatory part, then the optional part, and delay as much as possible the beginning of the action part, while of course executing it before the task deadline.

The scheduling algorithm used in the real-time control layer is based on the slack time server algorithm , as its characteristics are close to the REAKT requirements:

• The slack time server algorithm original goal is to maximise the CPU allocation to soft aperiodic tasks, which is more or less equivalent to maximise the amount of time available for optional processing. • The basic principle of the algorithm is to delay as much as possible the execution of periodic tasks while guaranteeing their deadlines. This is exactly what needs to be achieved in REAKT for action tasks. • The algorithm is an optimal one for the deadline-monotonic priority assignment, in the sense that no other algorithms can preserve more time for optional processing. • The slack time server algorithm is able to recover CPU time unused by stochastic real-time tasks, i.e. tasks which execute faster than their worst-case execution time.

The slack time server algorithm must be adapted to the REAKT task model for two main reasons:

• In the original algorithm, periodic and aperiodic tasks are considered independent. This is not the case in the REAKT task model, as dependencies exist between the mandatory, optional and action parts of a task. • Existence of sporadic tasks is not taken into account in the original algorithm.

The slack time server approach is based on the creation of a high-priority process, the Slack Server, which is in charge of computing the amount of optional time (slack time) available at a given time. When slack time is available and some optional tasks are ready to execute, the Slack Server pre-empts all real-time tasks to let the lower priority optional tasks run. The Slack Server must precisely know the amount of slack time which is available, in order to resume in time the pre-empted tasks. The problem is then to determine the algorithms which allow the Slack Server to compute the available slack time at each instant. Three issues must be considered:

• The first one describes the basic principles of the slack time server algorithm when considering only a periodic process set. • The second part is devoted to the links between the Slack Server and the second layer of scheduling. • Finally, approaches to scheduling sporadic tasks are considered.

THE SLACK TIME SERVER ALGORITHM

Principles

Consider a set of N periodic tasks, τ 1 , τ 2 …τ N . Each task τ i (C i , T i , D i , O i ) is characterised by its worst-case execution time C i , its period T i , its deadline D i and its initial release time O i . This model is compatible with the REAKT task model, assuming that:

• only the mandatory and action parts of a task are executed, with an execution time equal to the sum of the mandatory and action part worst-case execution times,

• sporadic tasks behave like periodic ones, with a period equal to their minimal interarrival time (worst-case behaviour).

The tasks are indexed in decreasing priority order, with priorities assigned to tasks according to the deadline-monotonic algorithm [Leung and Whitehead 82], i.e. the task with the shortest deadline receives the highest priority. One additional requirement is that all periods in the system must be integral multiples of some time unit, thus allowing to determine the hyperperiod H of the system (H is the least common multiple of all periods in the system). The hyperperiod is the minimal time after which the task activation pattern of the system repeats itself.

In order to compute the slack time, i.e. the maximum delay that can be inflicted at a given time to the whole task set without missing any deadlines, it is first necessary to compute the maximum delay that can be inflicted to a task τ i .

At time t, the total workload W i (t) at priority level i is given by:

The workload W i (t) gives the cumulative work that has arrived for priority 1 to i in the time interval [O, t].

Now consider the jth execution of τ i in the hyperperiod: the job release time R ij is given by O i +(j-1)T i and the job deadline D ij is R ij + D i . A necessary and sufficient schedulability condition (given in [Lehoczky et al. 89]) when using the DM priority assignment is:

which is equivalent to:

(1)

The computation of the maximum value for the above function does not need to be done on the continuous interval [R ij , D ij ], but only on the set of level-i scheduling points for the interval. A level-i scheduling point is defined as a time at which one or more tasks of priority higher than i are released. The set of scheduling points SP ij for the jth job of τ i is given by:

As the workload W i (t) is an monotonous increasing step function (with jump points given by the scheduling points), the difference t -W i (t) is a piecewise increasing function, which reaches a local maximum just before a jump point. Therefore, it is correct to search the maximum value only at the scheduling points. In fact, it is also necessary to consider the value of t -W i (t) at the deadline D ij , if it does not coincide with a scheduling point.

From the above considerations and equation (1), it comes that the maximum delay S ij that can be inflicted to the task τ i before the completion of its jth job is given by:

(2)

S ij is the maximum value which can be added to the workload W i (t) while still verifying equation (1). The series of S ij is non-decreasing for a given i if the task set is schedulable.

However, S ij is not the amount of optional processing that can be performed without missing the deadline of the jth execution of τ i : instead, S ij is the amount of work that can be performed at priority level lower than i. This includes of course optional processing, but also computation time of lower priority tasks and possibly idle time (if no optional work is ready when slack time is available). Therefore, we can introduce the level-i inactivity I i (t) as:

where Cj(t) is the time spent in task τ j , O(t) is the amount of optional processing already performed and I(t) is the time during which the processor has been idle.

The level-i slack time (i.e. the maximum amount of optional processing which can be performed while meeting task τ i deadlines) is thus given at any time by:

(3)

where φ i (t) gives the job number of τ i at instant t.

The global slack time at instant t S(t) is obtained by taking the minimum value of all level-i slack time at that instant. This is the maximum value by which all tasks in the task set can be delayed and still meet their deadlines.

(4)

The formulas given above show that the computations required to obtain the slack time are relatively complex. Fortunately, most of these computations can be performed offline: the series S ij , obtained by applying equation (2), involve only compilation time information, such as periods, deadlines and worst-case execution times. Therefore, a two-dimension table (the Slack Table) can be built off-line with these values, allowing fast access to the information at runtime. One drawback of the Slack Server approach is that it may request a large amount of memory: the size of the Slack Table can become quite large, as the first dimension is the number of jobs of τ 1 in the hyperperiod (assuming τ 1 has the shortest period of all tasks, and therefore the largest number of jobs in the hyperperiod), and the second dimension is the number of tasks in the task set. In REAKT, given the amount of memory requested by some AI components, the size of the Slack Table is not considered a major limitation.

Runtime behaviour

At runtime, the basic Slack Server algorithm is as follows:

1. compute the available slack time;

2. if slack time is available, pre-empt all periodic processes to let optional tasks run. After all the available slack time has elapsed, resume all periodic processes;

3. wait until more slack time is available and return to step 1.

The two main points to be discussed in the above algorithm are:

• the slack time computation;

• the time at which more slack time becomes available.

Equation (4) gives a formula to compute the slack time at instant t. In order to transform this formula into an algorithm, it is necessary to introduce a few counters:

• Idle will hold the cumulated idle time since the beginning of the hyperperiod;

• Optional will hold the cumulated optional processing time since the beginning of the hyperperiod; • C i will hold the cumulated time spent in task τ i since the beginning of the hyperperiod.

• n i will hold the number of completed executions of task τ i in the hyperperiod.

All counters are reset to 0 at the beginning of each hyperperiod. Each time the processor finishes an activity (idle, optional or periodic processing), the time spent in the activity is added to the corresponding counter. Each time a task τ i completes an execution, the counter n i is incremented.

The algorithm is the following:

LET inactivity = Idle + Optional min_slack_time = infinity // a large value

FOR i = Ntasks to 1 DO slack_time = slack_table[i][ni] -inactivity IF slack_time < min_slack_time THEN min_slack_time = slack_time ENDIF inactivity = inactivity + Ci ENDDO RETURN min_slack_time
This algorithm is a simple minimum value computation. The main point to note here is that the computation of the level-i inactivity I i (t) is performed incrementally by starting the minimum value computation from the lowest priority process, thus saving an iteration. The complexity of this algorithm is O(n).

The second point to study is the time at which more slack time becomes available. Once the currently available slack time is exhausted, more slack time will become available only after the blocking task (i.e. the task which was responsible for the minimum slack time value in the previous slack time computation) finishes.

Proof: Let b be the priority level responsible for the minimum slack time value S b in the last slack time computation. After exhaustion of this slack time, all level-i slack times S i are reduced by S b (after equation (3), since the optional processing time has been increased by S b ). In particular, available slack time at level b is 0. Available slack (which is the minimum of all level-i slack times) will remain 0 until the level-b slack time increases. Equation (3) shows that level-b slack time can increase only when S bj increases (as level-i inactivity can never decrease), i.e. when the job number of τ b changes. Therefore, it is necessary to recompute the slack time only after the completion of the current execution of τ b .

Stochastic tasks

As mentioned previously, the slack time server algorithm can be adapted to the recovery of unused CPU time after the completion of a stochastic task.

When stochastic execution occurs, the workload at priority level i W i (t) is overestimated, since it is based on the static worst-case execution time of the tasks. The real workload W' i (t) is given by: where C i (t) denotes the real cumulated execution time of task τ i (including expected computation times of task executions already started but no yet completed), WC i is the static worst-case execution time of task τ i and ΔC i (t) its cumulated saved time in completed executions. It is necessary to consider only completed executions, as there is no way to know how much will be saved in an on-going or future task execution.

The real maximum delay S' ij which can be inflicted to the jth job of task τ i is given by:

The above formula uses information (W' i (t) ) which can only be known at runtime, but it is clear that computing dynamically the series S' ij is not a practical solution. Therefore, an approximation of the formula is needed.

Let P ij be the first scheduling point in SP ij for which the maximum value S ij is obtained, and let the additional slack time Δ i (t) be:

Δ i (t) is a non-decreasing function, therefore:

From the definition of P ij :

From the above properties, it can be deduced that:

In order to consider S ij +Δ i (t) a valid approximation for S' ij (it will remain a valid approximation as long as it remains below the exact value of S' ij ), it is necessary to prove that the formula will be used only for values of t smaller than P ij The formula is used only for the jth execution of task τ i . Therefore, if it can be shown that the jth job of τ i completes before P ij , the approximation will be valid.

The proof is by contradiction: let t 1 be the time Δ i (t) was last computed (t 1 < P ij ), and suppose that:

The above condition states that the total workload (periodic workload plus computed delay) of the system at level i is higher than the available time, and thus that the jth job of τ i does not complete before P ij . Since t ≥ t 1 , Δ i (t) ≥ Δ i (t 1 ) . The above equations implies that: This is in contradiction with the definition of P ij : Equation (3) can now be updated with the above approximation:

(3')

The algorithm used to compute the available slack time must be modified according to the above formulas. It is first useful to define a new set of counters ΔC i which will hold the cumulated saved time since the beginning of the hyperperiod. The counters are reset at the beginning of each hyperperiod and the relevant counter is updated at the end of each task execution.

The algorithm for computing the slack time is now:

LET inactivity = Idle + Optional, min_slack_time = infinity FOR i = Ntasks to 1 DO slack_time = slack_table[i][ni] -inactivity IF slack_time < min_slack_time THEN min_slack_time = slack_time ENDIF min_slack_time = min_slack_time + ΔCi inactivity = inactivity + Ci ENDDO RETURN min_slack_time
As in the previous algorithm, it is important to note that the computation of the additional slack time Δ i (t) is done incrementally and does not require an additional iteration. The complexity of the algorithm remains O(n).

The property that additional slack time becomes available only when the blocking task completes is not true any more in case of stochastic executions: it is clear in this case that the early completion of a non-blocking task allows the system to recover some additional slack time. Therefore, the slack time computation algorithm must be executed after the completion of each process finishing early.

LINKS WITH THE SECOND SCHEDULING LAYER

Principles of interactions

In order to understand the interaction between the Slack Server and the REAKT components used in the second scheduling layer, it is useful to first list all relevant processes present in the REAKT Kernel, following decreasing priorities:

• The process with the highest priority is the Slack Server.

• Periodic and sporadic real-time tasks τ i , i = 1…N, sorted by decreasing priorities.

• The Intention Scheduler schedules the optional tasks (intentions) according to their deadline and importance. The complete scheduling algorithm used by the Intention Scheduler is beyond the scope of this paper. It is based on the earliest-deadline-first algorithm, and takes into account extensions proposed by [Chetto et al. 90]. • Optional parts of tasks are implemented as intentions.

All the activities which execute at a priority level below those of the real-time tasks can only run while these real-time tasks are suspended by the Slack Server. The basic system cycle is the following: • When slack is available, the Slack Server suspends itself and all real-time tasks;

• The Intention Scheduler is activated and can be in the following states:

♦ No intentions are ready to execute (all intentions are pending): Control is immediately returned to the Slack Server, which will resume the real-time tasks. This allows to preserve the slack time for later optional processing. The Slack Server needs to be activated when new information (e.g. events) able to trigger optional processing is sent to the Intention Scheduler, in order to compute the currently available slack time and let the optional processes run. ♦ Some new intentions are ready and need to be inserted in the schedule: This insertion is performed following the second-layer scheduling algorithm, using a function returning the amount of slack time available between now and the deadline of the intention. The specification of this function is given below. Once intentions have been scheduled, the algorithm can proceed to the following step. ♦ The schedule is not empty but the first intention is not active: This happens when the mandatory part associated to the intention has not been executed yet. In this case, the Intention Scheduler requests the Slack Server to resume only the specified mandatory task. Explicitly resuming real-time work from the second scheduling layer is the way of implementing at runtime the precedence constraints of the REAKT task model. ♦ The first intention in the schedule is active: it is executed while slack time is available. • When the slack time is exhausted, the Slack Server resumes all real-time tasks, thus pre-empting lower priority activities.

The interactions between the Slack Server and other REAKT components (mainly the Intention Scheduler) require two types of modifications in the slack time server algorithms:

• The basic loop of the Slack Server process must be modified, in order to react to messages sent by the Intention Scheduler. • The slack time computation algorithms must be extended to take into account the possibility of executing real-time tasks while slack time is available: this can occur when the Intention Scheduler requests the execution of the mandatory part of a ready intention.

Modified Slack Server loop

The Intention Scheduler can send three types of messages to the Slack Server:

• A notification that no more optional work is ready. The Slack Server should then resume all real-time tasks, in order to preserve available slack time for later use. This message can only be received while slack time is available, i.e. while the Slack Server is sleeping, waiting for the time-out to arrive. • A notification that new optional work is ready. This is the opposite message, and the Slack Server should then immediately recompute the available slack time and preempt the real-time tasks to start the optional work. This message can only be received while the Slack Server is not sleeping. Messages of this type should be generated only when the Intention Scheduler is idle and either an event is received or a pending intention becomes ready. • A request for the execution of the mandatory part of a given task. The Slack Server should then recompute the amount of slack time available, taking into account the early execution of the real-time task, and suspend itself for the newly computed duration. This message can only be received while slack time is available.

The above classification shows that steps 2 and 3 of the algorithm presented in section 5.2 should be modified. The modified algorithm is:

1. compute the available slack time;

2. if slack time is available, pre-empt all periodic processes to let optional tasks run. After all the available slack time has elapsed or a message is received from the Intention Scheduler, resume all periodic processes; 3. wait until either more slack time or optional work becomes available and return to step 1.

Early execution of mandatory parts

When the Slack Server receives a message from the Intention Scheduler requesting the execution of the mandatory part of a task, it must:

• resume the given task;

• compute the new available slack time, taking into account that the resumed task will execute first. Intuitively, the available slack time (i.e. the time available for executing both the resumed task and optional work) should increase, as the operation is basically an CPU exchange between optional processing and the real-time task. However, the exact amount of additional slack time depends on the priority of the resumed task.

Let r be the priority level of the resumed task. From a scheduling standpoint, executing τ r explicitly while slack time is available is equivalent to considering that the current execution of the task is optional processing and that the computation time of the mandatory part is null (this assumption remains valid while slack time is available). In turn, this implies that the workload for priority levels lower than or equal to r has been overestimated. The real workload W' i (t) is given by:

Therefore, the real level-i slack time S' i (t) is given by:

The above formula can be straightforwardly derived from the usual slack time formulas (equations ( 2) and ( 3)), as the additional term in the equation is constant (it has no effect on the maximum value computation involved in the computation of the series S ij ).

Therefore, the loop to compute the available slack time when a real-time task is explicitly resumed is simply divided in two parts: the first part computes the level-i slack for priority levels lower than or equal to r, adding C r to the slack time, the second for priority levels higher than r.

This algorithm has several interesting properties:

• When the resumed task has a priority lower than that of the blocking process, it is easy to show that no additional slack is available: the minimum value for the slack time, which is given by the level-b slack time (where b is the priority level of the blocking process), remains unchanged, as r > b. • When the resumed task has the highest priority, the additional slack time is equal to the task computation time: since all level-i slack time are incremented by the task computation time, it is clear that the minimum value is also incremented by this amount. • The property above can be usefully applied to real-time data acquisition tasks. These tasks usually have the highest priorities in the system. In addition, it is clear that they should not be suspended by the Slack Server, since all other tasks (including optional ones) use the data values they produce. Therefore, when a data acquisition task becomes ready while slack time is available, it should be immediately resumed. The above property can then be applied to compute the new available slack time, which should be the last computed slack time incremented by the data acquisition task computation time.

Computation of the available slack time in a given time interval

When scheduling optional work, the Intention Scheduler must take into account the computational requirements of real-time tasks: the available time to schedule an intention is equal to the amount of optional processing which can be performed between the intention release time (usually the current time) and the intention deadline.

The approach to compute this value is obviously related to the slack time computation: indeed, the available optional time O is given by: where D is the intention deadline.

However, it is difficult to apply the above formula directly: while it relatively easy to obtain the value of I i (t) at a given instant, by appropriately updating counters, the computation of this value in the future (which is the case for D) requires to know the number of times each task will complete and be released during the interval [t, D].

Computing dynamically this value would be too costly.

The solution used in REAKT is the following:

• compute the function giving the available slack at any time, assuming that optional processing is always ready. The values of this step function can be computed off-line and stored in an array for each time point in the hyperperiod. • the available slack time is then given by adding the currently available slack time to the difference between the value of the function at D and the current value of the function. By keeping appropriate time counters, it is possible to compute the available slack time until D without performing any iterations.

One aspect which is not considered in the current version, but on which more work is needed, is the existence of periodic intentions (optional parts of periodic tasks). Clearly, when scheduling an intention, the system should be aware of the workload that will be created by the forthcoming releases of periodic intentions. An off-line computation of a static periodic intention schedule, to be updated dynamically, could be an interesting approach.

SCHEDULING SPORADIC TASKS

The general problem of scheduling jointly periodic, sporadic and optional tasks is a difficult one, for which no general solution has been given in the literature: the underlying problem is to determine exactly the maximum amount of optional processing that can be performed, while the exact pattern of sporadic task arrival is by definition unknown. Therefore, there is either the risk to perform too much optional processing, and thus miss hard deadlines if many sporadic tasks are released, or to be too conservative and never perform any optional processing.

Two approaches have been studied for REAKT:

• The first one consists in applying a period transformation to the sporadic tasks, in order to transform them in stochastic periodic ones • The second one consists in computing the Slack Table assuming that the sporadic tasks behave like periodic ones, and then recover unused CPU time at runtime.

Period transformation

This approach consists in transforming a sporadic task τ i with deadline D i and computation time C i into an equivalent task with deadline 1 and computation time C i /D i . In this case, the task can be considered as a periodic task which is not always released, and the technique described for stochastic periodic tasks applies.

However, this approach has several drawbacks:

• It generates an overhead, as each task must be controlled by a timer in order to execute exactly for the planned amount of time.

• It may be difficult to apply when tasks use critical sections, as it is usually a good practice to avoid suspending tasks inside a critical section. • It increases the load of the system, as the tasks are executed at a higher priority than their original one (as their new deadline is smaller than before).

For all these reasons, it has been decided not to implement this solution.

Slack time recovery

In the second approach, the Slack Table is computed assuming the worst-case behaviour of the system, i.e. all sporadic tasks are considered periodic. The problem is therefore to recover slack time when the sporadic tasks are not executed. It is different from the problem of recovering unused time in stochastic task executions, as the sporadic tasks are not released are predefined instants, but rather can start at any time, as long as they respect their minimum interarrival time.

A sporadic task σ i of priority level i affects the value of S j (t), j ≥ i, in two ways:

• the computation of the workload W j (t),

• the computation of the set of scheduling points SP jk .

as both computations depend on the number of releases and the release times of σ i .

In order to avoid the complete slack time computation at runtime, it is necessary to find an approximation of the slack time which can make use of the values stored in the Slack Table . This in turn requires that the set of scheduling points used to compute each value in the table remains unchanged, as any modification would induce a change in the maximum value computation in equation ( 2).

The sets SP jk ,j ≥ i, remain unchanged when the release time of σ i is translated by T i , the minimal interarrival time for the task. Therefore, the approximation is to consider that slack time can be recovered as soon as the sporadic task has been delayed by T i , the amount of recoverable slack time being equal to the worst-case computation time C i .

The algorithm defined for dealing with stochastic task execution can then be used to compute the modified slack time.

This approach has been implemented, and experiments have shown that the approximation is very pessimistic, especially in the case of sporadic tasks having large interarrival times. For applications where time constraints are not too strict, dynamic incremental computation of the Slack Table could provide a better solution.

CONCLUSIONS

We have presented in this paper some elements of the REAKT Kernel, a software architecture well-adapted to the requirements of real-time knowledge-based applications. The paper focused on the scheduling techniques which allow the combination of predictable and unbounded computation while guaranteeing response times. The approach is based on the progressive reasoning paradigm, which allows a system to obtain a first solution quickly, and then refine it while time is available.

Activities in the REAKT Kernel are represented with the REAKT task model, which decomposes each activity into three parts, mandatory, optional and action. This model supports the easy integration of reflex behaviour and refinement activities, which is essential in the progressive reasoning approach. A scheduling algorithm based on the slack time server approach allows the system to maximise available time for optional processing while guaranteeing response times of the real-time tasks. This algorithm provides an answer to the complex problem of jointly scheduling periodic, sporadic and optional tasks in a real-time system.

A second scheduling layer, not described in this paper, is used to control the execution of optional processing. This layer can use both complex scheduling techniques, based on deadlines, importance and solution quality, and complex unbounded problem-solving algorithms. Interactions between the two layers is supported, and allows the second scheduling layer to modify, within certain limits, the behaviour of the real-time layer when optional processing requires it.

Validation of the REAKT Kernel has been performed through the development of several applications, including a large alarm management system in an oil raffinery. The development of this particular application has already shown that the REAKT functionalities are well-adapted to the requirements of a typical RTKBS such as a process control system.


Thierry A Mara
email: mara@univ-reunion.fr

Noura Fajraoui

Anis Younes

Frederick Delay

Inversion and uncertainty of highly parameterized models in a Bayesian framework by sampling the maximal conditional posterior distribution of parameters *

Keywords: Inverse modeling, Bayesian parameter estimation, Model parameter identification, Highly parameterized model, Heterogeneous transmissivity field, Karhunen-Loève expansion

We introduce the concept of Maximal Conditional Posterior Distribution (MCPD) to assess the uncertainty of model parameters in a Bayesian framework. Although, Markov Chains Monte Carlo (MCMC) methods are particularly suited for this task, they become challenging with highly parameterized nonlinear models. The MCPD represents the conditional probability distribution function of a given parameter knowing that the other parameters maximize the conditional posterior density function. Unlike MCMC which accepts or rejects solutions sampled in the parameter space, MCPD is calculated through several optimization processes. Model inversion using MCPD algorithm is particularly useful for highly parameterized problems because calculations are independent. Consequently, they can be evaluated simultaneously with a multi-core computer. In the present work, the MCPD approach is applied to invert a 2D stochastic groundwater flow problem where

Introduction

Models are tools on which environmental risk-assessment and decisionmaking strategies can rely, provided it is proved that the models are relevant to the problem under investigation. This relevance can be addressed by facing a model prediction to observation data knowing that the whole procedure also requires assigning model parameter values. Some parameters can be directly measured while some others ought to be indirectly estimated by comparing model predictions with observations. The present work addresses the issue of parameter identification for highly parameterized models. The notion of identification encompasses seeking the parameter values and assessing the uncertainty on parameters and on model predictions.

During the past two decades, the increasing power of computers was conducive to emphasize and promote the so-called Bayesian parameter estimation techniques. In essence, the Bayesian framework leads to the definition of the parameter joint posterior probability density function (pdf), for instance inferred by means of Markov Chain Monte Carlo (MCMC) samplings ( [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte carlo sampling methods using markov chains and their applications[END_REF][START_REF] Bayarri | A framework for validation of computer models[END_REF][START_REF] Higdon | Computer model calibration using high-dimensional output[END_REF]).

The notion of posterior pdf is associated with the fact that the parameter's pdf is conditioned both on plausible (prior) parameter values and on observation data. MCMC provides draws directly sampled from the posterior pdf which leads to exploration of the plausible areas in the parameter space.

The Bayesian estimation using MCMC has been subject to many developments and improvements during the last decade (e.g. [START_REF] Haario | An adaptive metropolis algorithm[END_REF][START_REF] Green | Delayed rejection in reversible jump metropolis-hastings[END_REF][START_REF] Haario | Dram: Efficient adaptive mcmc[END_REF][START_REF] Ter Braak | Differential evolution markov chain with snooker updater and fewer chains[END_REF] among others).

However, MCMC samplers remain computationally expensive because many draws are rejected by the statistical test embedded in the sampler. Furthermore, with MCMC, the parameters marginal posterior distributions cannot be investigated independently. Recently, several strategies have been proposed to increase MCMC efficiency (see [START_REF] Kennedy | Bayesian calibration of computer models[END_REF][START_REF] Christen | Mcmc using an approximation[END_REF][START_REF] Ma | An efficient bayesian inference approach to inverse problems based on an adaptive sparse grid collocation method[END_REF][START_REF] Cui | Bayesian calibration of a largescale geothermal reservoir model by a new adaptive delayed acceptance metropolis-hastings algorithm[END_REF][START_REF] Laloy | Efficient posterior exploration of a high-dimensional groundwater model from two-stage markov chain monte carlo simulation and polynomial chaos expansion[END_REF]).

In the present work we propose a new method, partly grounded in optimization techniques, to cope with the identification of model parameters.

The first step of this approach is to seek all the probable local optima of the joint posterior pdf of the whole set of parameters (including the maximum a posteriori estimate). Next, several maximizations of the conditional pdf are performed for different prescribed values of one selected parameter. The values assigned to this parameter are picked from a range around its optimal value(s). The value of the other parameters is investigated by maximizing the conditional pdf. This provides what we call the Maximal Conditional Posterior Distribution (MCPD) of the selected parameter. It actually corresponds to a discrete approximation of the pdf of a single parameter conditioned on data such that the conditional pdf is maximized.

The MCPD returns information about the model parameter values supported by the data and any correlations between parameters. The MCPD sample also allows uncertainty bounds to be assigned to the model predictions. The main advantage of the approach is that MCPD inferences for different parameters are independent and can be evaluated simultaneously by easily distributing the calculations over a multi-core computer (or several computers). This feature drastically decreases the computation time and makes the inversion of highly parameterized problems feasible.

The main topics addressed in the present paper are organized as follows. A short outline on inverse modeling within a Bayesian framework is proposed in Section 2, and then followed by the details on the MCPD sampling in Section 3. The first exercise testing the MCPD approach is proposed in Section 4 and addresses the ability of the sampler to retrieve a multimodal probability density function. The second test in Section 5 applies the MPCD approach to identify the Karhunen-Loève expansion ( [START_REF] Loève | Probability Theory[END_REF]) of a stochastic transmissivity field for a two-dimensional steady-state groundwater flow problem.

Bayesian inference

In inverse modeling, the parameter set (of size s) θ = {θ 1 , . . . , θ s } of a given model is estimated from a set of observation data d. In the following, we assume that the model does not suffer from misconceptions. The model is therefore supposed to be exact regarding the processes and the system that it mimics. However, observation data remain uncertain (random variables) making the model parameters to be also random and characterized by a joint probability density function p(θ). We denote by Ω i the probable prior uncertainty range of θ i . In a Bayesian framework, the parameter joint posterior pdf is defined by

p(θ|d) = p(d|θ)p(θ) p(d) (1)
where p(d) is a scaling factor called evidence, p(θ) is the prior density corresponding to a first guess on parameters before collecting the observations, while p(d|θ) is called the likelihood function measuring how well the model describes the data.

The parameter set that maximizes Eq. ( 1)

θ M AP = arg max θ p(θ|d) (2)
is called the maximum a posteriori estimate. It is the most probable parameter set given our knowledge about the system (i.e. the data d and the prior pdf of the parameters p(θ)) and it is sought by appropriate optimization algorithms (e.g., descent methods, evolutionary algorithms, etc...). Unfortunately, finding θ M AP does not allow to (fully) characterize the posterior uncertainty of the parameters (except for linear models, see [START_REF] Tarantola | Inverse Problem. Theory and Methods for Model Parameter Estimation[END_REF]). This uncertainty should be assessed by calculating the marginal posterior density for each parameter, defined as follows

p(θ i |d) = p(θ i , θ -i |d)dθ -i , ∀i = 1, . . . , s (3)
where θ -i represents the vector of parameters θ without θ i . The integral in

(3) can be approximated by a multidimensional quadrature method or by a sampling-based method such as the Markov Chain Monte Carlo (MCMC).

Nevertheless, the computational effort can be prohibitive and sometimes unaffordable for problems with a large number of parameters.

In the present work, we propose an optimization-based method in order to assess the parameter uncertainty for models post-conditioned on available observation data. For this purpose, we introduce the concept of maximal conditional posterior distribution. One could raise that relying on an optimization-based method will require solving many problems, as is classical with standard inversion techniques when obtaining a large set of solutions is contemplated. As shown hereafter, the maximal conditional posterior distribution has some specific features diminishing the calculation loads.

Maximal conditional posterior distribution

The concept

We define the maximal conditional posterior distribution (MCPD) of θ i as

P i (θ i ) = max θ -i (p(θ -i |d, θ i )) × p(θ i |d) (4)
P i (θ i ) is interpreted as the posterior probability function that maximizes the conditional posterior distribution p(θ -i |d, θ i ) and encompasses the MAP probability (i.e. P i (θ

M AP i ) = p(θ M AP |d))
. By using the Bayes theorem, one can write, max (p(θ

-i |d, θ i )) × p(θ i |d) = max (p(θ i |θ -i , d) × p(θ -i |d)).
Therefore, the MCPD in (4) can also be viewed as the distribution of the parameter θ i , knowing that the other parameters θ -i are at their optimal values. The MCPD of θ i is assessed in a discrete form by sampling Eq. (4).

A parameter θ i is frozen at a prescribed value and the other parameters within its prior uncertainty range Ω i (see Fig. 1). This gives,

θ * -i = arg max θ -i p(θ -i |d, θ i = θ * i ) (5)
P i (θ * i ) = p(θ * -i |d, θ * i ) × p(θ * i |d) = p(θ * |d) (6)
On the one hand, if θ i is globally identifiable, one expects that the further the parameter value is from θ M AP i , the more

P i (θ * i ) < P i (θ M AP i
). On the other hand, if θ i is not identifiable, varying its sampled values θ * i , will not change the joint posterior distribution in [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF]. Then, max (p(θ -i |d, θ i )) = p(θ M AP -i |d) and the MCPD of θ i will be equal to its marginal prior distribution (i.e. P i (θ i ) ∝ p(θ i ), see Eq. ( 4)). This is the case for instance of the unimodal multi-Gausian probability density function.

In the event of multimodality of p(θ|d), θ M AP is not unique or its search is hampered by the existence of local maxima. It is advisable to first acquire all the possible maxima of p(θ|d) by starting optimization processes from different locations in the parameter space (a multi-start procedure). Note that these prior searches do not ensure that all maxima will be identified.

When the local maxima are known, the MCPD calculation principle evoked above is generalized by sampling θ * i in all the subareas of the parameter space enclosing the local maxima. Finally, note that, by definition, the MCPD P i (θ i ) and the marginal posterior density p(θ i |d) are different. However, these densities would be exactly the same whenever max (p(θ

-i |d, θ i )) is constant ∀θ i ∈ Ω i .
It can be questioned on the relevance of the MCPD if the latter is not the marginal posterior density. We remind that the aim of the MCPD is to evaluate uncertainties on both the parameters and the model predictions given some observation data. The key point is less to know whether one can approach the marginal posterior density rather than providing uncertainties for valuable solutions to the inverse problem. As told above, by construction an MCPD samples the distribution of a parameter θ i knowing that the others are optimal. This new way brought by the MCPD to envision the uncertainty associated with a parameter seems to be an interesting alternative to more classical definitions. One can mention however that a Markov chain Monte Carlo sampler tries to draw all the probable solutions to the inverse problem while only some of these solutions are inferred by the MCPD sampling. In some situations, this can represent a drawback of the MCPD sampling.

MCPD assessment

For the sake of clarity, the algorithm inferring the MCPD assumes that the posterior density p(θ|d) has only one mode. The occurence of multimodal densities is just an adaptation of the procedure depicted below. The procedure starts by seeking the maximum a posteriori θ M AP (Eq. ( 2)). Then, the algorithm sketched by the flowchart in Fig. 2 is launched. Calculation of MCPD for each parameter takes place in two stages. The first stage identifies the relevant range within which the parameter θ i will be made to vary. For this purpose, we define a large sampling step, for θ i , e.g. ∆ × θ M AP i , with

∆ = 0.5. Then, θ i is successively set to θ ±k i = θ M AP i (1 ± k∆), k = 1, 2, 3, • • •
and the associated optimization in Eq. ( 5) is solved. Increments of k are stopped when:

1. the two current probability values p(θ +k |d) and p(θ -k |d) decrease below a prescribed threshold which may be expressed thus: p(θ ±k |d) < p(θ M AP |d)/100; or 2. the sampled values θ ±k i are beyond Ω i , i.e. the prescribed prior uncertainty range on θ i .

Usually, the first stopping criterion is reached before the second one, thus leading to a narrower sampling range of θ i compared to its prior uncertainty range Ω i . The second stage of the algorithm resamples the values of θ i to increase the refinement of the discrete MCPD. Before proceeding, the results from the first sampling step are re-ranked by increasing values of θ i . We

denote {θ k i , P i (θ k i i ) = p(θ k i |d), k i = 1, • • • , n i } the re-ranked first sample with θ k i i < θ k i +1 i .
We seek in this sample the interval where the difference

P i (θ k i +1 i ) -P i (θ k i i ) is maximal, i.e., we seek the index k m verifying k m = arg max k i P i (θ k i +1 i ) -P i (θ k i i ) (7)
Then, the optimization in Eq. ( 5) is solved for

θ * i = θ km+1 i + θ km i 2 and the pair (θ * i , P i (θ * i )) joins the set {θ k i , P i (θ k i i )}.
The latter is then sorted again by increasing value of θ i and the search of the index k m is resumed. This procedure makes that one resamples θ i and amends the set {θ k i , P i (θ k i i )} progressively (iteratively) at locations in the set where the MCPD is the most coarsely discretized. Usually, a few iterations on k m (i.e. N it = 10) are enough to obtain a good discrete depiction of the whole MCPD of θ i .

In the case of multimodal densities p(θ|d), the algorithm is adapted so that it repeats the first and second steps of sampling around each "optimal" θ opt value (see the example in Section 4).

At this stage, it is interesting to note that the optimization procedures to identify optimal vectors θ -i for prescribed values of θ i are completely independent of the optimizations identifying θ -j for prescribed values of θ j . Therefore, the MCPD samplings for each parameter of the investigated problem are independent and can be handled easily with parallel computing streams. There can be as many sessions as the number of parameters, which strongly reduces the total computation time. Obviously, other inverse techniques can be parallelized, for example, the Markov Chain Monte Carlo (see [START_REF] Braak | A markov chain monte carlo version of the genetic algorithm differential evolution: easy bayesian computing for real parameter spaces[END_REF], [START_REF] Vrugt | Accelerating markov chain monte carlo simulation by differential evolution with self-adaptive randomized subspace sampling[END_REF], [START_REF] Vrugt | Equifinality of formal (dream) and informal (glue) bayesian approaches in hydrologic modeling?[END_REF]). In this case, several chains are launched in parallel to explore the parameter space. The chains will generate independent subpopulations of solutions and exchange good individuals between the subpopulations to accelerate the convergence (the benefit of some emigration between subpopulations). The parallel calculations are no longer independent and a master computer is needed to analyze the independent subpopulations and generate new ones. In the end, the ease of parallelization brought by the MCPD inference is suited to the inversion of highly parameterized problems (see the example in Section 5).

Predictive density of an MCPD sample

Let us denote by n = s i=1 n i the total number of MCPD draws. A specific MCPD sample

(θ k , P i (θ k i )), k = 1, • • • , n
Given that, P i (θ k i )) = p(θ k |d), Eq. ( 8) is simply approximated by

p(d * |d) = n k=1 P i (θ k i )p(d * |θ k ) n k=1 P i (θ k i ) (9)
where p(d * |θ k ) is the likelihood function defined in Eq. ( 1) evaluated at d * and conditioned onto θ = θ k .

Inference of a multimodal distribution with MCPDs

The first case study deals with the ability of MCPD samplings to retrieve a three-modal probability density function with connected and disconnected modes. The test function to be retrieved is inspired from a case study proposed by Laloy and Vrugt who state that inferring a density function with disconnected modes by means of MCMC is very challenging ( [START_REF] Laloy | High-dimensional posterior exploration of hydrologic models using multiple-try dream (zs) and highperformance computing[END_REF]). In the present work, the density function also encloses modes that differ for each parameter. The density function to retrieve has 25 parameters and is the sum of three multiGaussian density functions expressed as

p(θ) = 1 6 N (µ 1 , 5C) + 2 6 N (µ 2 , 5I 25 ) + 3 6 N (µ 3 , 5I 25 ) (10)
where N (µ i , I 25 ) is the multiGaussian distribution of mean vector µ i . I 25 is the 25-dimensional identity matrix which indicates that the parameters (θ 1 , . . . , θ 25 ) are independent in the second and third Gaussian distributions in [START_REF] Christen | Mcmc using an approximation[END_REF]. C is a correlation matrix with null off-diagonal elements except for C 1,11 = C 11,1 = -0.5 and C 1,13 = C 13,1 = 0.8. These non-null terms impose, for the first Gaussian distribution in Eq. ( 10), a negative correlation between θ 1 and θ 11 and a strong positive correlation between θ 1 and θ 13 .

The three modes of each parameter are grouped in the vectors of means The first hurdle in this exercise is to identify the 3 local maxima of the density function. For all parameters, the maximal mean values are located in the third multiGaussian distribution with µ 3 as the vector of means. For this purpose we ran the optimization procedure twenty times, using initial solutions uniformly sampled from the parameter space. This calculation required about 6 500 evaluations of the density function (Eq. ( 10)). Subsequently, MCPD sampling was carried out in the vicinity of each optimum requiring an additional computation effort of about 3 000 runs for N it = 20 iterations of MCPDs refinements (see Section 3).

µ 1 = [-12, . . . , 12
In Fig. 3 we compare the maximal conditional posterior density of some parameters to their marginal posterior density. The MCPDs were computed numerically while the marginal posterior densities were assessed analytically by computing p(θ i ) = p(θ)dθ -i . Interestingly, when the modes are disconnected the two densities are equal (e.g. θ 1 ) whereas they can be very different when the modes are superimposed (e.g. θ 13 ). This is due to the fact that the marginal posterior pdf is an integral making that the superimposed modes are summed-up. Note that because the MCPDs are calculated by several optimization procedures in the vicinity of each mode, the superimposed modes produce overlapping MCPD curves.

It is noticeable that with only a few point estimates of the MCPD, the densities obtained are accurate enough (e.g. with only N it = 10 refinement iterations, not reported in Fig. 3). The "off-diagonal" scatterplots in Fig. 3 for a column i and a row j correspond to the pairwise MCPD draws (θ k i i , θ k i j ) and (θ

k j i , θ k j j
). The first one corresponds to the optimal sought value of the parameter θ j for sampled (prescribed) values of θ i and the second one to optimal values of θ i for sampled θ j . The first observation is that one can find several distinct values of a parameter θ i for the same draws of θ j . This is specific to the targeted pdf with overlapping modes. For example in Fig. 3, θ 13 has its three modes in (0, 13, 13), i.e. two overlapping modes, whereas modes of θ 1 at (-12, 1, 25) do not overlap. When seeking the MCPD of θ 13 , one identifies two local optima (θ 1 = 1, θ 13 = 13) and (θ 1 = 25, θ 13 = 13), i.e., several distinct values of θ 1 (around 1 and 25, respectively) are found for the same draws of θ 13 around the value 13.

The scatterplots also prove that the MCPD sampler is able to retrieve the correlation structure of the pdf. One can easily check that, in the vicinity of the first mode (for which θ 1 = 1, θ 11 = 11 and θ 13 = 13), a negative correlation between (θ 1 , θ 11 ) and a positive correlation between (θ 1 , θ 13 ) are observed. It can also be noted that no correlation is observed elsewhere, as assumed by the targeted density function (see Eq. ( 10)). Hence, the plots of θ 13 versus θ 11 draw orthogonal crosses (see scatterplot at row 3, column 2).

To conclude on this first exercise, we note that finding out the local optima is the most expensive stage. Identifying the multiple local optima is performed via a multi-start procedure launching searches from different initial locations in the parameter space. This procedures does not guarantee that all the local optima will be found. Our empirical experience shows however that the MCPD sampling (second stage) can identify some missed modes or improve the evaluation of some others. In this first stage, it is hard to state whether an optimization technique is a mark above the others. The use of gradient-based methods or evolutionary algorithms such as genetic algorithms or shuffled complex evolution methods are possible choices ( [START_REF] Duan | Optimal use of the sce-ua global optimization method for calibrating watershed models[END_REF]). For the second stage dealing with MCPDs sampling, the starting points in the parameter space are not far from the optimal locations. The calculations are fast and optimizations relying on gradient-based methods are recommended because they converge rapidly when started close to the solution. In the present example, the computation time was also strongly reduced by launching independent MCPDs samplings using 25 processors. This resulted in a cost of 260 computational time units (CTU) for the step identifying the local optima, whereas MCPDs sampling took only 120 CTU.

Inverting a stochastic field with MCPDs

We consider here an inverse problem with a two-dimensional random field as the unknown parameters.The two-dimensional random field of scalar values is denoted by Y (x, ω) with x the location in the Euclidean bounded space D and ω a realization index of the random field equivalent to a coordinate in the probability space Ω. Let us also consider a model response vector G(Y (x, ω)) and some observation data d. The inverse problem consists of finding an optimal estimate Ŷ M AP of the random field and its uncertainty range given the data. At first glance, the problem is ill-posed because there are an infinite number of unknowns, or at least, as many unknowns as the number of grid cells discretizing the domain D for numerically evaluating G(Y ). Hence, the first task is to reduce the dimensionality (regarding parameters) of the problem.

Karhunen-Loève expansion

Several authors have suggested using Karhunen-Loève (KL) transformation of the random field ( [START_REF] Laloy | Efficient posterior exploration of a high-dimensional groundwater model from two-stage markov chain monte carlo simulation and polynomial chaos expansion[END_REF][START_REF] Li | Efficient geostatistical inverse methods for structured and unstructured grids[END_REF][START_REF] Efendiev | Preconditioning markov chain monte carlo simulations using coarse-scale models[END_REF][START_REF] Marzouk | Dimensionality reduction and polynomial chaos acceleration of bayesian inference in inverse problems[END_REF]) to reduce the dimensionality of the problem. It is assumed that the random field obeys a second-order stationary Gaussian process, making that Y (x, ω) ∼ GP (µ Y , C Y (x 1 , x 2 )), where µ Y is the mean of the process, (x 1 , x 2 ) is a pair of different locations in the Euclidean space D, and C Y (x 1 , x 2 ) is the two-point covariance of the random process. This covariance is a scalar continuous and positive-definite function

corresponding to C Y (x 1 , x 2 ) = E [(Y (x 1 ) -µ Y ) (Y (x 2 ) -µ Y )]
, with E[ ] the mathematical expectation. Provided that Y (x, ω) is a real-valued random field with finite second moments, its KL expansion is

Y (x, ω) = µ Y + +∞ i=1 λ i ξ i (ω)ϕ i (x) (11)
where the douplets

(λ i , ϕ i (x)), ∀i ∈ N * are a set of eigenvalues (λ 1 ≥ λ 2 ≥ • • • > 0)
and the associated eigenfunctions of the covariance kernel C Y while {ξ i } ∞ i=1 are independent Gaussian random variables of zero mean and unit variance (ξ i ∼ N (0, 1)). The doublets eigenvalues-eigenfunctions are obtained by solving the Fredholm equation

D C Y (x 1 , x 2 ) ϕ i (x 2 )dx 2 = λ i ϕ i (x 1 ) (12)
When ranked by decreasing values, the eigenvalues tend more or less rapidly to zero, thus allowing truncation of the KL development to the K th order

Ŷ (x, ω) = µ Y + K i=1 λ i ξ i (ω)ϕ i (x) ( 13
)
where K is chosen so that K i=1 λ i ≥ (1ǫ) +∞ i=1 λ i . Given that the eigenfunctions {ϕ i (x)} ∞ i=1 are continuous and form a complete orthonormal system in L 2 (D), i.e. D ϕ i (x)ϕ j (x)dx = δ ij with δ ij the Kronecker delta function, Eq. ( 13) is therefore an expansion of Y (x, ω) onto an orthogonal basis. Provided the parameters µ Y and C Y (•) are fixed (guided by previous investigations), identifying Ŷ (x, ω) essentially involves seeking the plausible Gaussian random variables ξ * that explain the observation data d. The dimensionality of the inverse problem is therefore reduced to K parameters, i.e. the number of eigenvalues-eigenvectors chosen to expand Y (x, ω). The problem is also partly regularized since ξ is a vector of independent Gaussian variables, i.e. p(ξ) = N (0, I K ). Note also that one can include the mean µ Y in the sought parameters, as done in the sequel where one defines the parameter vector

θ = (µ Y , ξ) of dimension s = K + 1.
If we compare the above parameterization of a multiGaussian field with, for example a pilot point method, the reduction of dimensionality when solving the inverse problem is of the same order. Obtaining several plausible solutions in a Bayesian framework will finally require several optimizations procedures. On the one hand, a classical pilot point technique will duplicate the optimization of the pilot point values to modify different prior guesses on the parameter field. On the other hand, calculating MCPDs for the parameters of the Karhunen-Loève decomposition will also require multiple optimizations (see Section 3). Further studies handling carefully designed comparisons would be needed to identify advantages and drawbacks of both approaches.

Problem statement

Let us now illustrate the approach mixing a Karhunen-Loève decomposition and MCPD calculations. It must be first emphasized that the Karhunen-Loève decomposition is well known to depict accurately any second-order stationary random field ( [START_REF] Efendiev | Preconditioning markov chain monte carlo simulations using coarse-scale models[END_REF][START_REF] Marzouk | Dimensionality reduction and polynomial chaos acceleration of bayesian inference in inverse problems[END_REF]). Along this line, our study does not provide new insights on the KL decomposition. The main aim here is to assess the capability of MCPD samplings to correctly identify a large number of parameters that do not directly enter in the forward problem. As is the case with a lot of parameterization techniques, the optimized parameters serve as seeds to recompose the parameter field that enters in the calculation of the forward problem. In addition, the parameters of a Karhunen-Loève decomposition before inversion are independent because the eigenvectors of the decomposition are orthogonal. After inversion and post-conditioning, these parameters may become correlated. It appears interesting to see how MCPD samplings will cope with this modification and which kind of uncertainty can be derived for parameters that are the seeds of a parameterization technique.

Let us now illustrate the approach by considering the following steadystate two-dimensional flow problem in a heterogeneous medium

             ∇ • (T (x)∇h(x)) = 0, x = (x, y) ∈ D = [0, 1] × [0, 1] ∂h(x) ∂y y=0 = ∂h(x) ∂y y=1 = 0 T (x) ∂h(x) ∂x x=0 = 0.02; h(1, y) = 0 (

Posterior distribution

First, let us assume, as a prior guess, that the log-transmissivity field Y (x) has a mean of zero and an isotropic Exponential covariance C y (x 1 , x 2 ) = 10e -8|x 1 -x 2 | yielding an effective correlation length of 0.375. This choice of an erroneous covariance function compared with that of the reference problem is twofold. First, we want to show that it is still possible to obtain a good estimation of the random field even though the prior guess on its spatial structure is flawed. For instance, the correlation length of the guess is here half the correlation length of the reference field, and the prior variance is five times the variance of the reference field. One can also mention to justify this discrepancy that the model of covariance in practical applications is often more conjectured than really known. Second, the Fredholm equation in [START_REF] Cui | Bayesian calibration of a largescale geothermal reservoir model by a new adaptive delayed acceptance metropolis-hastings algorithm[END_REF] has well-known analytical solutions for an exponential covariance kernel (see [START_REF] Zhang | An efficient high-order perturbation approach for flow in random porous media via karhunen-loève and polynomial expansions[END_REF]) which simplify and accelerate the calculations of the eigenvalues-eigenfunctions. It is of course possible to choose another covariance kernel. In that case, one can make use of the calculation method proposed in [START_REF] Phoon | Implementation of karhunen-loeve expansion for simulation using a wavelet-galerkin scheme[END_REF] to estimate the eigenmodes. In the present study, the KL expansion is performed by keeping the first 103 modes representing 84% of the variance associated with the field Y (x). In the end, 104 parameters are sought (including the mean µ Y ), which is still a highly parameterized problem but far less than the initial problem that had a number of parameters equal to the number of meshes (10 4 ) discretizing the domain.

Given the way the 50 local data were obtained, the likelihood function is written

p(d|θ, σ m h , σ m Y ) ∝ exp - SS h (θ) 2(σ m h ) 2 - SS Y (θ) 2(σ m Y ) 2 (15)
where SS h and SS Y are the sum of squared errors on the hydraulic head and log-transmissivity respectively. In the present inversion exercise, we cannot include (σ h , σ Y ) in the parameters to be estimated. The optimization algorithm overfits the data because the number of unknowns (s = 104) is much greater than the number of data (i.e. 50). Hence, they are fixed to their exact presumed value.

Keeping in mind that in the KL expansion, ξ ∼ N (0, I K ) and assuming a uniform prior density function for µ Y allows to write (by using Eq. ( 15)) the parameter joint posterior density function as

p(θ|d, σ m h , σ m Y ) ∝ exp - SS h (θ) 2(σ m h ) 2 - SS Y (θ) 2(σ m Y ) 2 K i=1 exp - ξ 2 i 2 (16)
Maximizing the above density function amounts to minimize the following weighted sum of squares,

J(θ) = SS h (θ) (σ m h ) 2 + SS Y (θ) (σ m Y ) 2 + K i=1 ξ 2 i ( 17
)
The parameters ξ being normally distributed a priori, we did not prescribe any restrictive variation ranges. Regarding the parameter µ Y , we prescribed a very large variation range [-100, 100] making in practice µ Y free of any constraint. The type of objective function in Eq. ( 17) is rapidly minimized with gradient-based methods if the Jacobian matrix ∂h(x, θ) ∂θ is calculated accurately ( [START_REF] Kaczmaryk | Improving dual-porosity-medium approaches to account for karstic flow in a fractured limestone. application to the automatic inversion of hydraulic interference tests[END_REF][START_REF] Ackerer | Inversion of a set of well-test interferences in a fractured limestone aquifer by using an automatic multi-scale parameterization technique[END_REF]). A close look at the forward problem in [START_REF] Loève | Probability Theory[END_REF] shows that its differentiation with respect to any parameter θ parameter yields

             ∇ • (T (x, θ)∇h θ (x, θ)) = -∇ • (T θ (x, θ)∇h(x, θ)) , x ∈ D ∂h θ (x, θ) ∂y y=0 = ∂h θ (x, θ) ∂y y=1 = 0 T (x, θ) ∂h θ (x, θ) ∂x x=0 = -T θ (x, θ) ∂h(x, θ) ∂x x=0 ; h θ (1, y) = 0 (18)
where the notations h θ (x, θ), T θ (x, θ) stand for the derivatives ∂h(x, θ) ∂θ , ∂T (x, θ) ∂θ , respectively. Provided the forward problem has been solved and one can analytically differentiate the terms T θ (x, θ), the calculation of ( 18)

is very similar to that of the forward problem [START_REF] Loève | Probability Theory[END_REF]. It can be handled with the same code, the same discretization, resulting in comparable accuracy for evaluating both h(x, θ) and h θ (x, θ).

values. Though the patches of high and low values are correctly located, the kriged map differs from the reference field (Fig. 4(a)) and the MAP solution (Fig. 4(c)). For example, the effective correlation length seems overestimated and the overall map is smoother than the reference field. As shown later, the sampled variogram of the kriged map differs from that of MCPD solutions and confirms the preceeding visual appraisals. Stated differently, the local transmissivity values do not conceal all the information and it makes sense inverting a flow problem conditioned on head data to retrieve the logtransmissivity field. It can also be noted that the inverse problems solved to calculate the MCPDs overfit the data, the number of parameters being larger than the conditioning data. The consequence is that the terms enclosed in the objective function Eq. ( 17) become lower than the prior measurement errors on heads and transmissivity data.

The uncertainty on the predicted log-transmissivity (or hydraulic head) can be estimated with Eq. ( 9). The uncertainty (the 95% confidence interval size) associated with the log-transmissivity field is shown in Fig. 4(d). The crosses mark the (co)locations of heads and transmissivities data used for the inversion. We can note that the uncertainty is not constant over the domain and is lower at the measurement locations where the 95% confidence interval size is about 0.4 (i.e. ∼ ±1.96σ m Y ). This makes sense because the model prediction was anticipated to be more accurate at the measurement locations. Moreover, most of the domain is assigned an uncertainty range less than 0.8 (∆ log(T ) < 0.8) which is quite narrow except in the areas that contain no data (see North and South-West boundaries of the domain).

It is interesting to check whether the parameter uncertainty inferred via

MCPDs yields a set of log-transmissivity fields including the reference (true) field. Fig. 4(e) maps the areas (in black) where the local value Y (x) of the reference field is outside the 95% probable uncertainty range calculated by MCPDs (the width of these local uncertainty ranges are given in Fig. 4(d)).

As expected, the areas where the uncertainty ranges do not include the local values of the reference field are those where no prior information was collected to solve the inverse problem. In the present case however, the uncertainty ranges are relatively narrow (Fig. 4(e)) and the non-matching areas represent 30% of the total domain. They reach 13% for the 99% predicted uncertainty range. Interestignly, the uncertainty bounds encompass the true field in the vicinity of observed data where the uncertainty range is smaller (∆ log(T ) < 0.8). This is the consequence of a well-conditioned problem. Observation data are rather evenly spread over the whole domain and very few sub-areas are not documented by at least a transmissivity value giving the order of magnitude of the parameters. Finally, it is worth noting that increasing the variance of the postulated Exponential covariance kernel (e.g. σ 2 Y = 20) provides slightly wider uncertainty bounds that finally encompass the true field (not shown). This artificial increase of the σ 2 Y value to obtain increased uncertainties on the transmissivity field is not a consequence of flawed MCPD samplings. The handled problem is here highly regularized by the Karhunen-Loève decomposition which controls the freedom of action onto the tuned parameters (i.e. the eigenvalues of the decomposition).

The MCPD estimates of the first KL coefficients are reported in Fig. 5 ("diagonal" plots). We notice that the density of µ Y has a bell-shaped curve centered on zero, though a prior uniform density was postulated. Fig. 5 also reports, in its "off-diagonal" plots, on the various draws of parameters θ to build MCPDs. A plot in column i and row j corresponds to the optimal values of the parameter θ j for sampled (prescribed) values of θ i and optimal values of θ i for sampled θ j . These plots for the pair (µ Y , ξ 1 ) yield a unique straight line of non-null slope indicating a strong negative correlation between these two parameters. Notably, µ Y and ξ 1 are also closely correlated with ξ 6 (row #7, columns #1-2 of Fig. 5) meaning that these parameters have dependent effects on the state variable calculated by the flow model. We remind that the prior guess on the parameters ξ assumed that they were independent and of Gausssian distribution with zero mean and unit variance. The posterior estimates show that some parameters can be highly correlated and that their distribution may have changed. For example, the parameters from ξ 2 to ξ 4 show bell-shaped MCPDs, but are clearly non-centered on zero and have a variance less than one (width of the bell-shaped curve far less than 6).

The overall quality of the inverse solutions is now discussed. First, the MAP solution reported in Fig. 4(c) has a spatial structure close to that of the reference solution. Note that the reference solution has an effective correlation length of about 0.5 which is smaller than expected from the covariance model used to generate the random field. This occurs when the size of the domain is close to the correlation length. The sampled covariance (variogram) of the MAP solution is Gaussian with a correlation length of about 0.45, close to that of the reference field (see Fig. 6(a)). The MAP solution, however, slightly underestimates the variance, probably because the KL expansion considers modes representing only 84% of the variance on Y (see above). However, the MAP solution is much closer from the reference field than a kriged map of transmissivity data. This kriged map has a variogram that largely overestimates the correlation length and underestimates the variance of the transmissivity field. As told above the kriged field is smoother than the reference and MAP solutions. We note that some inverse solutions calculated by taking parameters from MCPDs sampling can perfectly fit the covariance of the reference field (Fig. 6(a)). These solutions are good but not the best with respect to the posterior pdf.

Note that MCPDs sample give optimal inverse solutions for a prescribed value of one parameter. Using the sample, we rebuilt the associated Y fields and calculated their sampled (experimental) covariances. Then, these covariances were fitted with a Gaussian model in the form of C y (x 1 , x 2 ) = σ 2 y e -η(x 1 -x 2 ) 2 , the variance of C Y obviously being σ 2 y and the effective correlation length 3/η. Fig. 6(b), (c), and (d) show these variances and correlation lengths. In general, the MCPD solutions slightly underestimate the variance of Y with mean values on the order of 1.9 for a reference at 2.1.

The mean correlation length of MCPD solutions establishes at 0.43 for a reference at 0.47. Though the prior guess on the spatial structure of the inverse solutions greatly differed from the reference, the eigenvectors using the KL expansion clearly record the effective spatial structure (covariance) from which data are extracted. The re-composition of these eigenvectors by means of MCPDs provides a large set of probable solutions that are only very slightly biased, at least for the forward problem handled here.

Conclusions

In this work, we introduced the concept of maximal conditional posterior distribution for the identification of model parameters. The MCPD of a parameter θ i can be viewed as the distribution of the parameter θ i knowing that the other parameters are at their optimal values (in the sense that the conditional posterior pdf is maximized). It makes sense calculating a MCPD because the uncertainty drawn from the distribution of θ i refers to valuable solutions to the inverse problem. All the parameters θ -i are optimal except θ i which is allowed to wander in the parameter space. In essence, the calculations of two MCPDs for parameters θ i and θ j are independent.

They can be performed very easily on parallel computational sessions and one can extend the number of parallel sessions up to the total number of parameters involved in the problem under investigation (thus, drastically reducing the computational effort). This feature associated with the fact that the MCPD draws are good inverse solutions post-conditioned on observation data, renders the inversion in a Bayesian context affordable, even for highly parameterized problems.

As stated above, a MCPD allows a single parameter θ i to vary when the others are optimal. A key feature to avoid the single parameter (and the others) to be far from valuable solutions is to pre-identify the so called maximum a posteriory (MAP). This MAP is a solution where all the parameters are at their optimal value, but this solution may correspond to several local maxima when the underlying distribution of parameters is multimodal.

Hence, the success of MCPD sampling depends on the ability of the preliminary optimization procedure to retrieve all the probable local optima. This can be achieved by multi-start optimizations. The synthetic test case performed in this study showed that the MCPD calculations retrieve fairly well multimodal distributions.

The MCPD technique was also faced to the problem of retrieving a parameter random field on the basis of information gathering both measurements of the parameters and state variables from a spatially distributed problem (here, steady-state groundwater flow). The parameterization allowing a strong decrease in the dimensionality (in parameters) of the problem was grounded in a Karhunen-Loève decomposition of the parameter field. This technique is well-known to be accurate regarding the depiction of a stationary random field. Hence, the MCPD calculations were really tested on their ability to provide valuable inverse solutions and the associated uncertainties for pa-rameters that do not directly enter in the calculation of the forward problem (the parameters allow the re-composition of a transmissivity field then used in the forward groundwater flow problem). The whole procedure allowed us to approximate the random field with only s = 104 parameters: the mean of the random field and the 103 first eigenmodes of the Karhunen-Loève decomposition. Then, the MCPDs of these parameters were assessed quickly with the proposed algorithm taking advantage of their computation in parallel.

It was shown that the MCPD-KL association was able to accurately retrieve a reference field even when starting from a flawed initial estimate of the spatial structure (covariance) of the field. MCPDs can also render parameter distributions that strongly differ from their prior guess especially regarding their mean and variance. Set k i = 2, ∆ = 0.5, it = 1, P i (θ 1 i ) = p(θ M AP |d)

If k i even, θ k i i = θ M AP i (1 -k i 2 ∆) Else, θ k i i = θ M AP i (1 + (k i -1) 2 ∆)
Compute θ k i -i from Eq.( 5) Set P i (θ k i i ) = p(θ k i |d)

k i = k i + 1 i = i + 1 max(P i (θ k i i ), P i (θ k i -1 i )) < p(θ M AP

θ

  -i are optimized to find (according to the Bayesian definition) the maximal probability of these parameters. Changing the prescribed value of θ i allows scanning the distribution of θ i . In practice, the sampled values of θ i (denoted below θ * i ) are picked around the MAP estimate θ M AP i (estimated beforehand)

  corresponds to a point location on the hypersurface described by the parameter joint probability density p(θ|d). The density function of the model predictions conditioned on observed data d allows assigning uncertainty bounds to the prediction of a new observation d * . If we assume that d * is independent of d conditional on θ (i.e. p(d * |d, θ) = p(d * |θ), see [19]), we can write p(d * |d) = p(d * |θ)p(θ|d)dθ

  ], µ 2 = [1, . . . , 25] and µ 3 = [25, . . . , 1]. Thus, θ 13 has two modes located at θ 13 = 13. To compute the MCPDs, and specifically to seek optimal vector θ * -i knowing the value of θ * i , we use the MATLAB optimization toolbox, especially the fminunc.m function based on the socalled trust-region method. This function minimizeslog(p(θ)) and the convergence velocity is accelerated by providing the values of the derivatives oflog(p(θ)) with respect to the θ components.

  Fig. 4(d) respectively.

  Figure captions

Figure 1 :

 1 Figure 1:

θ km i +θ km+1 i 2 Find

 2 Figure 2: 36

Figure 4

 4 Figure 4:

Figure 5 Figure 6 :

 56 Figure 5:

Acknowledgement

The authors are grateful to the French National Research Agency who funded this work (project RESAIN n • ANR-12-BS06-0010-02).

Results and discussion

With 104 coefficients kept in the KL expansion, MCDPs were calculated on an eight-core computer, each core calculating 13 MCPDs. On average, each core ran the model (forward flow problem + calculation of the Jacobian matrix) about 6 500 times (CTU) for a total over the eight cores of 54 560 runs. 105 spatially distributed problems were solved, one being the calculation of the head variable (Eq. ( 14)), and 104 evaluating sensitivities to parameters θ -i for a prescribed value of θ i (Eq. ( 18)). As a comparison, in [START_REF] Laloy | Efficient posterior exploration of a high-dimensional groundwater model from two-stage markov chain monte carlo simulation and polynomial chaos expansion[END_REF], the authors handled a similar problem by means of MCMC and following the approach proposed in [START_REF] Cui | Bayesian calibration of a largescale geothermal reservoir model by a new adaptive delayed acceptance metropolis-hastings algorithm[END_REF]. To lower the computational load, a surrogate model was used in a two-stage approach in which a new candidate must first pass the surrogate likelihood successfully before undergoing the statistical test of the original model. Three dependent chains were launched in parallel.

The chains started to converge after 4 000 calls of the forward (i.e. original) model for a total of 10 000 calls (CTU).

The reference log-transmissivity field is reported in Fig. 4(a). The logtransmissivity field stemming from the KL expansion using the maximum a posteriori estimate of the parameter vector θ M AP = (µ Y , ξ) M AP is reported in Fig. 4(c). Both fields in 4(a) and 4(c) closely resemble each other, despite the fact that the covariance kernel of the KL expansion differs from that of the reference field (see above). Note also that we deliberately chose an erroneous covariance kernel for the KL expansion. However, the sampled covariances of the reference field and of the MAP estimates are very similar, implying that initial errors on the type of covariance and on the correlation length can be amended by inverting the KL coefficients and mixing eigenfunctions to recompose a relevant random field.