Julien Erny

Josette Pastor

Henri Prade

A Similarity and Fuzzy Logic-Based Approach to Cerebral Categorisation

This work proposes a formal modelling of categorisation processes attempting at simulating the way information is categorised by neural populations in the human brain. The formalism mainly relies on a similarity-based approach to categorisation. It involves weighted rules that use inference and fusion techniques borrowed from fuzzy logic. The approach is illustrated by a simulation of the McGurck effect where the combination of contradictory auditory and visual stimuli creates an auditory perceptive illusion.

INTRODUCTION

The understanding and the prediction of the clinical outcomes of focal or degenerative cerebral lesions, as well as the assessment of rehabilitation procedures, necessitates knowing the cerebral substratum of cognitive or sensorimotor functions. This is addressed by functional neuroimaging studies. They have shown that sensorimotor or cognitive functions are the offspring of the activity of large-scale networks of anatomically connected cerebral regions [START_REF] Bressler | Large-scale cortical networks and cognition[END_REF]. However, no one-to-one correspondence between activated networks and functions can be found [START_REF] Démonet | A pet study of cognitive strategies in normal subjects during language tasks. influence of phonetic ambiguity and sequence processing on phoneme monitoring[END_REF].

To understand the intricate relations between these networks and the cognitive functions, we need a model dealing explicitly with symbolic information while being strongly rooted in our knowledge of the neurobiological processes.Two complementary approaches address the problem. On the one hand, connexionist models focus on neural representation, at different scales. Biologically-plausible modelling deals with individual neurons, achieving a good level of accuracy but losing the symbolic aspect of information [START_REF] Tiesinga | Computational model of carbachol-induced delta, theta, and gamma oscillations in the hippocampus[END_REF][START_REF] Wang | Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model[END_REF]. On a coarser scale, simpler neuron models have revealed links between the neuronal substratum and simple cognitive functions [START_REF] Arbib | Synthetic brain imaging: grasping, mirror neurons and imitations[END_REF][START_REF] Grossberg | Thalomocortical dynamics of the McCollough effect: boundary-surface alignment through perpeptual learning[END_REF][START_REF] Lumer | Neural dynamics in a model of the thalamocortical system. part 1. layers, loops, and the emergence of fast synchronous rhythms[END_REF][START_REF] Lumer | Neural dynamics in a model of the thalamocortical system. part 2. the role of neural synchrony tested through perturbations of spike timing[END_REF]. On the other hand, models from symbolic AI tend to reproduce cognitive processes without being interested in biological plausibility [START_REF] Anderson | A spreading activation theory of memory[END_REF][START_REF] Ghahramani | Modular decomposition in visuomotor learning[END_REF][START_REF] Minsky | A framework for representing knowledge[END_REF][START_REF] Minsky | The Society of the Mind[END_REF][START_REF] Taatgen | Modeling parallelization and flexibility improvements in skill acquisition: from dual tasks to complex dynamic skills[END_REF][START_REF] Taatgen | Why do children learn to say "broke"? a model of learning past tense without feedback[END_REF]. However none of these approaches seems to answer the question of how the activation of large-scale networks follows the brain's information processing. This paper proposes an intermediary layer of interpretation, which echoes some other works [START_REF] Arbib | Synthetic brain imaging: grasping, mirror neurons and imitations[END_REF][START_REF] Balkenius | Nonmonotonic inferences in neural networks[END_REF][START_REF] Labatut | Cerebral modeling and dynamic bayesian networks[END_REF]. However these works may lack a certain biological plausibility especially in the representation of categorisation processes, problem that will be addressed hereafter.

In the present paper, we first discuss the previous models and describe the requirements that a new model should address. Then we present a new approach that can fulfill some of these prerequisites. We then apply this model to phonemic discrimination. We conclude with a discussion of the results and some perspectives.

EXISTING MODELLING APPROACHES

Representing information processing in large-scale cerebral networks necessitates departing from traditional connectionist methods, which tend to represent neural populations as neurons [START_REF] Alexander | Do cortical and basal ganglionic motor area use "motor programs" to control movement?[END_REF].At this level, symbolic information processing is specific to the neuronal population that implements it and can be assumed to emerge from the numerical processing in the population's individual neurons. The causal connectivity approach [START_REF] Pastor | Information processing in large-scale cerebral networks: The causal connectivity approach[END_REF] integrates this symbolic processing with more numerical aspects linked to the activation levels of cerebral regions in large-scale networks. This approach inspired probabilistic models [START_REF] Labatut | Cerebral modeling and dynamic bayesian networks[END_REF] and has influenced our own work.

In causal connectivity, the neuronal populations are seen as information processors while the axon bundles linking these populations are information transmitters. The information itself is represented as two-dimensioned data: i) a numeric component, called magnitude, stands for the overall activation of the neuronal population that processed this piece of information (thus allowing comparisons with neuroimaging data), and ii) a symbolic component, called type, that qualifies the pattern of the firing neurons in the population. Information is then spread in the network of populations using a dynamic Bayesian network formalism. The treatment taking place in a node can be seen as categorisation since the type of incoming information is compared to archetypes stored in the node. These archetypes represent the configuration of the synaptic weights inside the population and are modified dynamically every time information passes through, thus implementing an unsupervised, on-line learning mechanism. However this model is lacking a primordial feature of cerebral categorisation: similarity [START_REF] Sloutsky | The role of similarity in the development of categorization[END_REF]. It is known that close-by neurons will activate and fire for similar stimuli, suggesting a similarity code between stimuli based on the proximity between neurons responding to them. In the model mentioned above [START_REF] Labatut | Cerebral modeling and dynamic bayesian networks[END_REF], there is no such a similarity and the archetypes are disconnected from each other. As a consequence, the model cannot interpolate between archetypes, losing the ability to compute a fine-tuned response to an ambiguous signal.

The information flowing through the brain is noisy and often uncertain. This is a major issue that has to be dealt with. Moreover, the proposed framework should have good potentials for learning on an unsupervised and online basis. Amidst classical categorisation systems in AI [START_REF] Michie | Machine Learning, Neural and Statistical Classification[END_REF], none meets exactly these requirements. However, case-based reasoning [START_REF] Aamodt | Case-based reasoning: Foundational issues, methodological variations and system approaches[END_REF][START_REF] Bergmann | On the role of abstraction in case-based reasoning[END_REF][START_REF] Bouchon-Meunier | Fuzzy sets and possibility theory in approximate and plausible reasoning[END_REF], and fuzzy rules-based reasoning [START_REF] Bouchon-Meunier | Fuzzy sets and possibility theory in approximate and plausible reasoning[END_REF] provide ideas that can be a starting point for a specific formalism.

Our goal is to build a model dealing explicitely with cerebral information, while being as close as possible to neurophysiological reality, and flexible enough to incorporate new knowledge coming from neurosciences. However, learning will not be addressed in this paper.

A MODEL OF A NODE IN A LARGE-SCALE CEREBRAL NETWORK

A node implements the information processing and propagation by a functionally coherent neuronal population. These mechanisms are based on a functional decomposition of this population.

Principle

Decomposition of a node

A node receives information through several inputs and emits one output. For each input, we consider a receiving subpopulation and, symmetrically, we consider another emitting subpopulation. Inside each subpopulation we consider patterns representing, i) in the case of the receiving populations, the neurons responding to a specific configuration of the input, and ii) in the case of the emitting population, the neurons firing to propagate the output. The output domain is discrete since the internal regions of the brain does not handle continuity (every continuous incoming signal, like audio waves for example, is "broken" into categories by primary sensory cortex). Receiving and emitting patterns are linked through inference rules using the following principle : an emitting pattern is linked to one receiving pattern in each receiving population. Hence, there is one inference rule per output category (i.e. emitting pattern). Although a given receiving pattern may be used in different rules, its contribution to the activation of the emitting pattern is rule-dependent. Therefore, what is referred to in the rule is not the receiving pattern per se (i.e. the neuron's bodies) but a functional image of this pattern as it is used in the rule (i.e. the activated neurons : bodies and axons). From now on, "receiving pattern" will be used to designate those functional images (see Fig. 1).

Figure 1:

Physical view: the neuronal population is partitioned in two receiving subpopulations and an emitting one. Each receiving population contains two receiving patterns connected to the emitting patterns (some of the the connections are not shown for readibility's sake). The functional view shows how two functional images of the same pattern can be used in two different rules.

Patterns

The patterns are formally fuzzy sets. In the case of emitting patterns, there exists one fuzzy set corresponding to each element of the output domain, that element being the core of the set. This element is called referent and characterizes the emitting pattern. The other members of the set will be referents of other emitting patterns and they belong to the set with a degree that reflects the similarity between these other patterns and the considered pattern. As for the receiving patterns, their cores contain the elements of the corresponding input domain that are perfectly recognised while the other members are the closeby and less recognised elements. The next section gives a formalism for these notions and presents the key equations describing it.

Inference Principle

Similarity of incoming information is computed against the receiving patterns. The more they are compatible, the more these patterns are activated. When the activity of a pattern exceeds a threshold, the pattern fires. When all the receiving patterns of an emitting pattern fire, this pattern is activated. Taking into account the similarity of the activated emitting patterns, a decision policy is then applied to decide which pattern fires the final output (see Fig. 2).

Formalism

Let X be a node with n inputs. Information is a couple (M, T ) where M is the magnitude and T , the type. M ∈ R + and T is a fuzzy set defined on D T , a discrete set. Let (M in i , T in i ) be the incoming information on input i and (M out , T out ) the output information. To keep the notations manageable, we call Di and Dout the domains of T in i and T out . Let BR = {Ri} i∈ [1,p] be the inference rule base. We call P R j i the fuzzy set standing for the receiving pattern in input i for the rule Rj, and P E j the fuzzy set standing for the emitting pattern corresponding to the conclusion of Rj . We call na j i (resp. sa j i ) the level of activation (resp. the threshold of activation) of the receiving pattern P R j i . We call dec j i the firing level, i.e max(0, na j isa j i ). Finally, we call dec j the level of activation of the emitting pattern P E j .

Activation during reception

For each time step, for each rule Rj and for each input i, the compatibility

α j i of T in i with P R j i is computed as max C∈D i min(T in i (C), P R j i (C))
, which estimates the consistency of the two fuzzy sets. Hence the activation level at time t is computed from the previous activation level by increasing it when α j i is sufficiently large and the pattern has not been already fired recently:

na j i (t) = relax (1)
X .na j i (t -1) +

k.fact(α j i , Mi(t)).frefract(dec

j i (t -1))
where,

• fact is increasing from [0, 1] × [0, 1] on [0, 1], and prevents information with a too low magnitude from passing through. M is a normalized magnitude (∈ [0, 1]). fact should be a threshold function, like a sigmoïd for instance.

•

f ref ract is decreasing from [0, 1] on [0, 1], with f ref rac (0) = 1 et f ref rac (1) = 0.
It stands for the refractory period of the population, preventing it to activate again when it has been activated strongly before. A linear function (f ref rac (dec) = 1dec) or again a sigmoïd can be used.

• relax

X ∈ [0, 1] is called the relaxation parameter and characterises the temporal integration of the pattern. If it equals 0, there is no temporal integration.

Activation during emission

We have to update dec j for each rule Rj. We assume that the combination of the emitting patterns obtained at the previous step is linear with respect to each compononent, hence,

dec j (t) = relax (2) X .dec j (t -1) + n i=1 pi.dec j i (t -1) (2)
where the pi are parameters that have to be tuned in order to produce the more accurate simulation, and where relax

X is the relaxation parameter for the emitting subpopulation, defined as above.

When the dec j are known, we still have to combine and decide which pattern will fire. If there is only one emitting pattern activated (dec k ≥ 0), then the task is done and T out = P E k . However, when several rules R k 1 , • • • , R kq are triggered, there is a conflict that has to be solved. But first we will cumulate the activations, meaning that every time a pattern is mentioned as close to another, they will share their activations. If rule(C) is the rule from which C is the referent, dec rule(C) is the total activation of the corresponding emitting pattern (i.e. direct activation summed with similarity activations):

dec rule(C) (t) = i∈{k 1 ,•••,kq } P E i (C).dec i (t) (3)
When all the activation levels are computed, we just choose the most activated pattern as the firing one.

Output magnitude

The magnitude representation is handled as causal connectivity as in [START_REF] Labatut | Cerebral modeling and dynamic bayesian networks[END_REF], by a Bayesian-like propagation mechanism. The output magnitude depends on the one hand on the nature of the node (its function) and on the other hand on the level of compatibility of the incoming information. Formally,

M out (t) = f (1)
X (decmax(t)).f

(2)

X (M in (t -1), u) (4)
where M in (t -1) is the average over all the incoming magnitudes at t -1; u is a random variable standing for the noise, and decmax is the activation of the firing emitting pattern. The function f (1) is increasing, defined from [0, 1] on [0, 1]. Since its purpose is to generate a low output magnitude for something badly recognised, a threshold function would be a good choice. f (2) is a real function that may be non-linear.

Example

Let us consider a very simple example to make these equations somehow clearer. The focus will be on the inference mechanism letting aside the temporal aspect. To do so, we will describe the behaviour over one time step t0, in a completely relaxed node (there are no activation carried over from previous steps) and with all activation thresholds reduced to 0 (i.e ∀i, j, na j i (t0 -1) = dec j i (t0 -1) = 0). Consider a node with two attributes 1 and 2, with domains D1 and D2, receiving respectively inputs T in 1 and T in 2 , with weights p1 and p2. The magnitudes for each input are set to 1 and the parameter k is also set to 1. There are three "if...then" rules in the base :

(P R 1 1 , p1) ⊗ (P R 1 2 , p2) -→ P E 1 (5) (P R 2 1 , p1) ⊗ (P R 2 2 , p2) -→ P E 2 (6) (P R 3 1 , p1) ⊗ (P R 3 2 , p2) -→ P E 3 (7)
Where ⊗ stands for conjunction. It then works in three steps :

• First, for each receiving pattern P R j i , reception activations na j i are computed. In this simplified case, it is just a matching rate between the pattern and the input (e.g. [START_REF] Bouchon-Meunier | Fuzzy sets and possibility theory in approximate and plausible reasoning[END_REF]), t in i being the input at t0: na j i (t0) = max

C∈D i min(T in i (C), P R j i (C)) (8)
• Then, the direct activation dec j of the emitting pattern P E j is computed :

dec j (t0) = i∈{1,2} pi.na j i (t0) (9)
• Finally, the direct activations of the emitting patterns are combined using similarity and according to equation (3). Fuzzy sets P E 1 , P E 2 and P E 3 are shown in Fig. 3. In the figure, c1 belonging to P E 1 (whose referent is c2) at a level h1 means that the similarity between c1 and c2 is h1. It also means that the activation of the pattern whose referent is c1 will contribute to the activation of P E 1 with h1 as a discounting factor. The equation (3) gives us: dec rule(c 1 ) (t0) = dec 2 +h2.dec 1 ; dec rule(c 2 ) (t0) = dec 1 + h2.dec 2 + h1.dec 3 ; dec rule(c 3 ) (t0) = dec 3 + h1.dec 1 , where dec rule(c i ) is the total activation of the pattern of referent ci. The most activated pattern (i.e. for which dec rule(c i ) is maximum) is then chosen as the final output.

APPLICATION TO PHONEMIC DISCRIMINATION

As an illustration we have built with a unique node a very simple phonemic discriminator. This has permitted to reproduce an auditory illusion, the McGurck Effect.

Model

To define a phoneme, we use the articulatory theory of language that characterises a phoneme by the way we produce it. Table 1 represents The node is built with already learned categories. Assume there are four receiving patterns : P R bil = {(bil, 1.0), (lab, 0.6)}, P R dea = {(lab, 0.6), (dea, 1.0), (alv, 0.6)}, P R vel = {(pal, 0.6), (vel, 1.0)} and P R plo = {(plo, 1.0), (plv, 0.6)}. There are three categories (i.e. emitting patterns) already implemented in the node : [b], [d], [g] and the corresponding rules are shown in expression [START_REF] Bressler | Large-scale cortical networks and cognition[END_REF]. The parameters are : k = 0.6, and relax (1) = 0.98 (for eq.1), and relax (2) = 0.15 (for eq.2).

McGurck Effect

Figure 5 shows the principle of this illusion discovered in 1976 by McGurck [START_REF] Mcgurck | Hearing lips and seeing voices[END_REF]. The visual [ga] is a person articulating a [ga] while the sound played is a [ba]. Using some recent insights about its neurophysiological nature, we implemented it using the same node but in a modulated way, i.e where the thresholds triggering the category [g] are lowered when a visual [g] is presented along with an auditory [b]. Thanks to the combination of activations in the emitting patterns, the response comes as a [d] (see Fig. 6). Because of the lowering of the activation threshold, the emitting pattern [g] gets activated, and the pattern [b] is activated by the stimulus itself. If we look at Table 1, we see that between these two sounds, stands the [d]. Hence, the pattern [d] is close to the other two and gets activated when they are, In this simple experiment, we have shown that this model of a node can discriminate language sounds, and even reproduce an auditory illusion that was impossible to reproduce in a simple manner with previous models based on causal connectivity. However, we are not claiming that this experiment is an accurate description of how the real McGurck effect is produced by the brain. At least, it shows a demonstration of how we can use the model to simulate observable phenomenons with a proper choice of parameters.

CONCLUSION AND FUTURE WORK

Cerebral modelling is particularly important in the medical field, to better understand the possible outcomes of a lesion or a degenerative disease. Causal connectivity [START_REF] Labatut | Cerebral modeling and dynamic bayesian networks[END_REF], although promising, is still lacking a certain biological plausibility. In this work we tried to fill this gap a little, and we managed, doing that, to reproduce a phenomenon that was not modelled before by the other models of this family.

Our approach found its roots in empirical knowledge on neurophysiology and we scaled it on an "intermediary level" : the neuronal populations. This interpretation finds an interesting echo and a beginning of formal justification in the work of Balkenius and Gardenförs [START_REF] Balkenius | Neural mechanisms for self-organization of emergent schemata, dynamical schema processing and semantic constraints satisfaction[END_REF][START_REF] Balkenius | Nonmonotonic inferences in neural networks[END_REF], who proved that it was possible to produce non-monotonic logical inferences using activation patterns of neuronal networks. The inferences could then be encoded by weighted formulas as in possibilistic logic [START_REF] Benferhat | Representing default rules in possibilistic logic[END_REF]. This calls for a comparison with the similaritybased inference mechanism as a future line of research.

Regarding the future work, dynamic learning mechanisms are being developed to build and tune the rules using only incoming information. In terms of machine learning, it is an online and unsupervised learning. We plan to test these learning abilities in two ways. They will be compared to other algorithms on machine-learning benchmarks. We will also design a simple categorisation experimental task and will compare results of the task performance by humans and by the model. Furthermore, using the work of Balkenius and Gardenförs as a starting point, we would like to formally prove that the neuronal subtratum can support high level logic when one reasons in terms of patterns. This would contribute to fill the gap between connectionism and symbolic logic. Lastly, it has been shown in biology that the action of different neurotransmitters modulates the behaviour of neuronal populations (e.g. threshold lowering as in McGurck effect, see 4.2). By modelling specifically this action, we would be able to modulate our system in a systematic manner, while being even closer to biological reality [START_REF] Loubinoux | A single dose of the serotonin agonist paroxetine enhances motor output: double-blind, placebo-controlled, fMRI study in healthy subject[END_REF].

  et al. (Eds.) IOS Press, 2006 © 2006 The authors. All rights reserved.

Figure 2 :

 2 Figure 2: Inference principle

Figure 3 :

 3 Figure 3: Emitting patterns

RFigure 4

 4 Figure 4 shows the response of the node to a stimulus [b]. The stimulus is presented during 50 ms, every second, to allow neuronal processing to take place. We see clearly that only the pattern [b] fires, the other two staying unactive all the time. Now what happens when we change the threshold?

Figure 4 :

 4 Figure 4: Activations of the emitting patterns, in response to an auditory stimulation of type [b]

Figure 5 :

 5 Figure 5: The simultaneous presentation of a 'visual' [ga] and of an 'auditory' [ba] leads to the auditory perception of a [da]

Figure 6 :

 6 Figure 6: Activations of the emitting patterns in response to an auditory stimulation of type [b], modulated by a visual stimulation of type [g].

Table 1

 1

	suggest and we only focus on [b], [d] and
	[g] (i.e. [b] "is close" to [d] and "less close" to [g] and [d] "is close"
	to [g]).

Table 1 :

 1 Consonants in the French phonologic system

	locus ⇒ mode ⇓	labial dental velar
	plosive	p	t	k
	voiceless
	plosive	b	d	g
	nasal	m	n

J. Erny et al. / A Similarity and Fuzzy Logic-Based Approach to Cerebral Categorisation


Akylas Stratigakos

Simon Camal

Andrea Michiorri

Georges Kariniotakis

Prescriptive Trees for Value-oriented Forecasting and Optimization: Applications on Storage Scheduling and Market Clearing

Keywords: Data-driven optimization, decision trees, electricity market, prescriptive analytics, value-oriented forecasting

Decision-making in the presence of contextual information is a ubiquitous problem in modern power systems. The typical data-decisions pipeline comprises several forecasting and optimization components deployed in sequence. However, the loss function employed during learning is only a proxy for task-specific costs (e.g. scheduling, trading). This work describes a data-driven alternative to improve prescriptive performance in conditional stochastic optimization problems based on nonparametric machine learning. Specifically, we propose prescriptive trees that minimize task-specific costs during learning, embedded with a scenario reduction procedure to reduce computations, and then derive a weighted Sample Average Approximation of the original problem. We present experimental results for two problems: storage scheduling for price arbitrage and networkconstrained stochastic market clearing, respectively associated with electricity price and load forecasting. The empirical results show significant improvements on deterministic and stochastic lookahead policies, with the relative difference being more pronounced when training data are limited.

I. INTRODUCTION

The ongoing transition towards the smart grid era is associated with an influx of data from various sources. Making informed decisions under uncertainty, entails leveraging such contextual information (also known as explanatory data or features) within the data-decisions pipeline. The conventional modeling approach comprises a sequential process, where forecasting models first estimate uncertain parameters (e.g. load, market prices), which are subsequently input into an optimization problem. However, learning is performed by minimizing a proxy loss function, without considering the impact of forecasts on downstream costs. Assessing this impact and designing algorithms that maximize forecast value, rather than accuracy, constitute key challenges in energy forecasting [START_REF] Hong | Energy forecasting: A review and outlook[END_REF] for the coming years.

Conditional stochastic optimization departs from the classical approach to study decision-making under uncertainty in Parts of this research were carried in the frame of the Smart4RES project (No. 864337), supported by the Horizon 2020 Framework Program. the presence of contextual information. The data-driven framework put forward in [START_REF] Hannah | Nonparametric density estimation for stochastic optimization with an observable state variable[END_REF] and extended in [START_REF] Bertsimas | From predictive to prescriptive analytics[END_REF] integrates machine learning with optimization by employing local learning (nonparametric) algorithms to infer a weighted Sample Average Approximation (SAA) of stochastic problems conditioned on contextual information. A related stream of research proposes learning under alternative loss functions to derive, possibly biased, forecasts that explicitly minimize downstream costs (or equivalently maximize forecast value) [START_REF] Elmachtoub | Smart "predict, then optimize[END_REF], termed cost-/value-oriented forecasts. Relevant applications on wind [START_REF] Li | Toward cost-oriented forecasting of wind power generation[END_REF] and load [START_REF] Zhang | Cost-oriented load forecasting[END_REF] forecasting rely on a two-step approach that involves first inferring a convex loss, then training the model, which, however, does not directly leverage the optimization component. In [START_REF] Morales | Prescribing net demand for electricity market clearing[END_REF] authors prescribe value-oriented demand forecasts to clear a day-ahead market, by considering forecast impact on balancing costs during learning. A generic endto-end learning approach to train probabilistic forecasting models to minimize task-specific costs for problems with a smooth objective is described in [START_REF] Donti | Task-based end-to-end model learning in stochastic optimization[END_REF], with applications in storage and grid scheduling. The trading cost for a risk-aware renewable producer is used as an alternative loss to forecast market quantities in [START_REF] Stratigakos | A value-oriented price forecasting approach to optimize trading of renewable generation[END_REF]. Lastly, empirical risk minimization (ERM) can be applied to directly forecast decisions, see e.g. data-driven solutions to the newsvendor problem [START_REF] Ban | The big data newsvendor: Practical insights from machine learning[END_REF], with relevant applications on renewable trading found in [START_REF] Munoz | Feature-driven improvement of renewable energy forecasting and trading[END_REF], [START_REF] Carriere | An integrated approach for valueoriented energy forecasting and data-driven decision-making application to renewable energy trading[END_REF], respectively involving linear decision rules and neural networks. A main drawback, however, is that the ERM approach fails to guarantee the out-of-sample feasibility of decisions. Considering a discrete set of actions circumvents this issue; see e.g. mode-based control of storage posed as a supervised classification problem [START_REF] Henri | A supervised machine learning approach to control energy storage devices[END_REF].

The work presented here integrates the framework established in [START_REF] Bertsimas | From predictive to prescriptive analytics[END_REF] with learning under alternative loss functions, describing a decision tree algorithm that directly minimizes task-specific costs during learning and outputs decisions rather than predictions (prescriptive trees), effectively learning a policy from data. We extend our previous work [START_REF] Stratigakos | Prescriptive Trees for Integrated Forecasting and Optimization Applied in Trading of Renewable Energy[END_REF] by considering a multi-temporal setting, describing an embedded scenario reduction procedure, and explicitly modeling recourse actions when generating prescriptions. To illustrate the proposed solution, we provide a comprehensive evaluation against deterministic and stochastic lookahead solutions derived under the standard modeling approach in two applications: price arbitrage with storage and network-constrained market clearing. Further, we provide empirical evidence on the effect of training sample size and regularization hyperparameters on the out-ofsample optimization performance. Overall, the main contribution of this paper is to propose and validate an effective data-driven alternative to the standard "forecast, then optimize" modeling approach, offering a generalized perspective that jointly examines value-oriented forecasts and decisions, directly leverages the underlying optimization problem, and ensures out-of-sample feasibility. The proposed approach is generic and readily applicable to convex conditional stochastic optimization problems.

The rest of the paper is organized as follows. Section II details the prescriptive trees methodology. Section III describes the applications and presents the results. Finally, we conclude and provide directions for future research in Section IV.

II. METHODOLOGY A. Prescriptive Analytics Problem

We consider optimization problems that involve uncertain parameters of interest ξ ∈ Ξ ⊆ R d ξ (e.g. load), associated with features γ ∈ Γ ⊆ R dγ (e.g. weather conditions). Our focus is on solving the prescriptive analytics problem

v = min z∈Z E P [c(z; ξ)|γ = γ] = min z∈Z E ξ∼Pγ [c(z; ξ)], (1)
where v is the objective value, z ∈ R dz is the decision vector, Z is a convex set of feasible solutions, c(•) is a cost function, γ = γ is a feature observation, P is the joint distribution of γ and ξ, and P γ is the marginal distribution of ξ conditioned on γ. Typically, P, P γ are unknown, thus we want to infer a solution from a training data set {(ξ i , γ i )} n i=1 of n observations, by learning a policy z that varies as a function of γ. In this work, ξ i typically represents a sample path realization of a stochastic process, e.g. one day of load observations. Further, we assume access to an oracle that efficiently generates an SAA of (1) with off-the-shelf solvers.

The standard modeling approach, termed "forecast, then optimize" (FO), suggests first deploying forecasting models to infer an approximation P γ of the marginal conditional distribution of ξ (or its expectation), and then solving

ẑF O = arg min z∈Z E ξ∼ Pγ [c(z; ξ)], (2)
which now defines a conventional stochastic optimization problem, that can be tackled with standard methods. However, as the forecasting and optimization components are separate, the decision-maker cannot effectively hedge against miscalibrated forecasting models. Integrating machine learning and optimization offers a viable alternative in this case. Following [START_REF] Bertsimas | From predictive to prescriptive analytics[END_REF], we consider a weighted SAA of (1) as

ẑP P (γ) = arg min z∈Z n i=1 ω n,i (γ)c(z; ξ i ), (3)
termed predictive prescription, where ω n,i (γ) ≥ 0 denotes weights learned from a local learning algorithm, e.g. nearest neighbors and decision trees, satisfying n i=1 ω n,i (γ) = 1. Recall that for the standard SAA ω n,i (γ) = 1 n ; thus we optimize over the empirical distribution of ξ without leveraging any contextual information. We extend this framework by proposing a methodology to derive ω n,i (γ) by directly minimizing expected costs [START_REF] Hong | Energy forecasting: A review and outlook[END_REF]. Following [START_REF] Mundru | Predictive and prescriptive methods in operations research and machine learning: an optimization approach[END_REF], we formally define the problem of searching over functions f : Γ → Ξ that improve prescriptive performance as

min f ∈ F, z f (γ i ) ∈ Z i∈[n] c(z f (γ i ); ξ i ) (4a) s.t. z f (γ i ) = arg min z∈Z i∈[n] ω f n,i (γ i )c(z; ξ i ) ∀i ∈ [n], (4b)
where [n] := {1, . . . , n}. In the following, we focus exclusively on decision tree ensembles, which generally outperform other local learning algorithms.

B. Prescriptive Trees

The proposed tree algorithm follows the popular classification and regression trees (CART) [START_REF] Breiman | Classification and regression trees[END_REF] method, which recursively applies locally optimal splits, leading to a set of L leaves, R 1:L . A node split defines a disjoint partition of feature space at feature j ∈ d γ and point s: R

1 (j, s) = {i ∈ [n]|γ ij < s} and R 2 (j, s) = {i ∈ [n]|γ ij ≥ s}.
Following [START_REF] Breiman | Classification and regression trees[END_REF], the locally optimal split that minimizes task-specific loss is derived from

min j,s   min z1∈Z γi∈R1(j,s) c(z 1 ; ξ i ) + min z2∈Z γi∈R2(j,s) c(z 2 ; ξ i )   .
(5) The inner minimization problems reflect the SAA solution of each partition, where z 1 , z 2 are locally constant decisions. Decision trees are generally prone to overfitting, which hinders their capacity. An effective way of mitigating the adverse effects of overfitting, i.e., high variance, is via randomizationbased ensemble methods such as Random Forests [START_REF] Breiman | Random forests[END_REF] and Ex-traTrees [START_REF] Geurts | Extremely randomized trees[END_REF], which can be readily adapted to our framework to grow ensembles of prescriptive trees (prescriptive forests).

Starting at the root node, we recursively partition the feature space until either a stopping criterion is met or performance cannot be further improved. Typical stopping criteria include the maximum tree depth ∆ max , the minimum number of observations n min per leaf, and a predefined threshold for cost reduction. At each tree node, we apply perturbation by selecting a subset of features and candidate split points (and possibly bootstrapping samples), which is an integral part of ensemble methods. The ExtraTree algorithm suggests randomly selecting K features and a random candidate split point per feature (without bootstrapping), while in the Random Forest algorithm all candidate splits per selected feature are evaluated (with bootstrapping). In this work, we opt for the ExtraTree algorithm due to its reduced computations (fewer splits are evaluated per node). After selecting K candidate splits, the aggregated cost [START_REF] Li | Toward cost-oriented forecasting of wind power generation[END_REF] for each one is estimated and compared with the cost at its parent node, updating the tree structure accordingly. Algorithm 1 details the subroutine of splitting a single tree node with a random split criterion.

Algorithm 1 NodeSplit

Input: Data {γ i , ξ i } n0 i=1 , current partition R 0 , current depth ∆, hyperparameters {n min , K, ∆ max } Output: Local split (j, s) 1: Estimate cost v 0 = min z∈Z γi∈R0 c(z; ξ i ) and prescription z 0 = arg min z∈Z γi∈R0 c(z; ξ i ) at current node R 0 2: if ∆ < ∆ max and n 0 ≥ 2n min then 3: Select K random features { γ 1 ,. . . , γ K } 4: Select random split point s κ ∀ κ ∈ K from U [γ min κ , γ max κ ] such that |R 1 |, |R 2 | ≥ n min 5:
Solve ( 5) ∀(κ, s κ ) {Optional: apply ReducedSAA} 6:

Find tuple

(j * , s * ) such that v(s * ) = min κ=1,••• ,K v(s κ ) 7: if v(s * ) < v 0 then 8: ∆ ←-∆ + 1 9:
return tuple (j * , s * ) {Create new split} 10:

else 11:
return nothing {Node becomes leaf with output z 0 } 12:

end if 13: end if Similarly to CART trees, for query γ the corresponding weights ω n,i (γ) are derived as:

ω n,i (γ) = I[γ i ∈ R(γ)] |R(γ)| , (6)
where R(γ) denotes the target leaf that γ belongs to, | • | denotes the set cardinality, and

I[γ i ∈ R(γ)
] is an indicator function to assess whether γ i falls into R(γ). For an ensemble of B trees the derived weights are

ω n,i (γ) = 1 B B b=1 I[γ i ∈ R b (γ)] |R b (γ)| . (7)
Single trees are fully compiled, i.e., leaves provide prescriptions without requiring further computations during implementation. For prescriptive forests, taking the average decision might lead to infeasibility, thus an additional step of solving (3) for a query γ is required.

Finally, it is straightforward to derive value-oriented forecasts as ξ = n i=1 ω n,i (γ)ξ i . The key difference is that the internal tree structure is now dictated by local splits that minimize the expected cost of a constrained problem, rather than forecast error, as implemented in the CART method.

C. Embedded Scenario Reduction

Problem ( 5) is of a discrete nature and generally intractable for all but the simplest problems, which motivates the application of the ExtraTree algorithm, as it requires fewer evaluations per node. We propose to further reduce computations by accelerating individual split evaluations by means of applying scenario reduction (the terms scenarios and sample paths are used interchangeably). First, we select an upper threshold for the number of scenarios n max , such that n max ≥ n min , which guarantees that scenario reduction will be applied in both partitions. Next, we apply a scenario reduction method, solve the SAA (5), and finally scale the expected cost back to aggregated cost. The process is detailed in Algorithm 2. Naturally, as we move down the tree structure, the average number of observations per node decreases. Therefore, scenario reduction mostly affects splits near the top of the tree, and is applied progressively less frequently moving downwards. We refer the reader to [START_REF] Dvorkin | Comparison of scenario reduction techniques for the stochastic unit commitment[END_REF], and references therein, for an overview of the literature on scenario reduction and a comparison of different methods to a stochastic unit commitment problem. For this work, we employ k-medoids clustering to select representative samples.

III. APPLICATIONS

To illustrate the efficacy of the proposed method, we compare it to the traditional modeling approach in two applications. The first considers scheduling a storage device to perform price arbitrage in a day-ahead (DA) market. Next, we consider stochastic market clearing with network constraints and uncertain demand. Thus, the two applications are associated with electricity price and load forecasting, respectively. To facilitate reproducibility of the results, we provide the code for the experiments in [20].

A. Forecasting Models

As discussed, the standard modeling approach requires first generating forecasts of uncertain parameters. Depending on the optimization module, these forecasts range from point forecasts, to probabilistic forecasts and trajectories (scenarios). As the purpose of this work is not to provide a comprehensive evaluation of forecasting models, but rather to showcase the efficacy of an alternative approach, we select established benchmarks. In the following, electricity prices and load forecasts are generated using a Quantile Regression Forests (QRF) model, which is a generalization of Random Forests. The QRF is robust to overfitting, requires little tuning, and can generate both point and probabilistic forecasts. Samplebased approximations of multi-period stochastic optimization problems require modeling temporal correlations between successive time periods; to this end, a Gaussian copula function is employed, given the conditional marginal predictive densities obtained from the QRF, following the procedure detailed in [START_REF] Pinson | From probabilistic forecasts to statistical scenarios of short-term wind power production[END_REF]. For the sake of completion, comparisons with naive benchmarks are also provided. For reference, Fig. 1 illustrates different types of load forecasts for the DA horizon.

B. Storage Scheduling and Price Forecasting

1) Problem description: The first application concerns scheduling a generic storage device, similar to [START_REF] Donti | Task-based end-to-end model learning in stochastic optimization[END_REF]. The operator of a grid-scale storage device decides the charging z in and discharging z out actions, which consequently induce a specific state of charge z soc , for each period t of the forecast horizon T . The operator wants to maximize profits via price arbitrage, while also considering battery degradation costs, and penalizing excessive deviations from the starting state. Here, we assume that degradation costs and excessive deviations (henceforth both referred to as penalties) are modeled as quadratic regularization terms. The problem, which depends on estimating the uncertain DA prices π da , is formulated as where η in , η out denote the charging and discharging efficiency, z 0 the initial state of charge, z (resp. z) the lower (upper) limit of decision variables, and [T ] := {1, . . . , T }. For simplicity, the state of charge must be at the starting level at the end of the day. Recourse actions are not considered. The following approaches (policies) are compared:

min z in t , z out t , z soc t E π da t∈[T ] π da t (z in t -z out t ) + γ z soc -z 0 2 2 + z out 2 2 + z in 2 2 (8a) s.t. z soc ≤ z soc t ≤ z soc ∀t ∈ [T ], (8b)
z in ≤ z in t ≤ z in ∀t ∈ [T ], (8c) z out ≤ z out t ≤ z out ∀t ∈ [T ], (8d)
z soc t+1 = z soc t - 1 η out z out t + η in z in t ∀t ∈ [T -1], (8e)
z soc 1 = z 0 , z soc T - 1 η out z out T + η in z in T = z 0 , (8f)
• FO-Deterministic (FO-Det): Generate point price forecasts, then solve a deterministic optimization problem, which is the standard modeling approach. • Prescriptive Forest (PF): Generate predictive prescriptions with the proposed prescriptive forest algorithm. This is equivalent to solving a deterministic optimization problem with value-oriented forecasts. • SAA: Determine an SAA of (8) given sample paths of {π DA } n i=1 , which minimizes in-sample cost. The SAA serves as a naive benchmark, since it ignores contextual information.

2) Experiment setup: This step employs a total of 3 years of market data from France, obtained from the ENTSO-E Transparency platform. Features γ include lagged prices, predicted net load at the system level, and categorical variables for month, weekday, and hour of the day. We consider a batchmode learning setting and vary the size of the training sample n from 6 months to 2 years to assess the sensitivity of the out-of-sample performance. We also examine the impact of regularization parameters γ, (for simplicity assumed equal). Table I presents the complete list of indicative parameter values used. For all cases, evaluation is performed on the last year of data. To train the PF, the hyperparameters are set as {B = 50, K = d x , n min = 10}, with individual trees being maximally grown, and without scenario reduction.

3) Results: Fig. 2 summarizes the out-of-sample prescriptive performance of all of the considered cases. Overall, the PF consistently outperforms the other benchmarks, resulting in the highest profit, while also setting the efficient frontier in all cases (values towards the top and left are better). With respect to the standard FO-Det, for a given value of γ, , the PF generally incurs similar penalties combined with significantly higher profits. Regarding sensitivity to sample size, in general, larger samples translate into improved optimization performance, which is expected as the learning component improves. Regularization penalties remain fairly stable throughout; however, attained profit steadily increases with the sample size, with results converging for 2 years of training data. The PF and the SAA exhibit less sensitivity to sample size compared to the FO-Det approach, which fails to outperform the naive SAA solution for the smallest sample examined (6 months). Thus the SAA poses an attractive alternative for the case of limited training data. For larger sample size, both the PF and the FO-Det outperform the SAA, with results being more pronounced for larger values of γ, . Finally, Fig. 3 examines forecast accuracy under a standard loss function (from FO-Det) and an alternative loss function (from PF), highlighting that, in this case study, improved accuracy does not translate into improved optimization performance.

C. Stochastic Market Clearing and Load Forecasting

1) Problem description: Next, we examine market clearing mechanisms in conjunction with the different modeling approaches. We consider clearing a forward (day-ahead or DA) market and a real-time (RT) market with network constraints and uncertain demand. In the DA horizon, the system operator sets the dispatch levels of controllable generators in order to satisfy the uncertain demand, while respecting a set of technical and network constraints. When uncertainty is realized during the RT operation, corrective dispatch actions are required to maintain the demand-supply equilibrium. A DC linear approximation of network constraints is applied, as is standard practice. For simplicity, the unit commitment problem (start-up/shut-down actions) is not modeled. Note that in this subsection we are primarily interested in improving the out-ofsample prescriptive performance for a given market clearing mechanism, rather than examining in detail its economic properties. The following approaches are compared:

• Deterministic market clearing under FO (FO-Det): Generate point estimates of uncertain load, then solve a deterministic problem to clear the network-constrained DA market. This defines a multi-period DC-OPF problem, and serves as the base solution.

• Stochastic market clearing under FO (FO-Stoch): Generate probabilistic load forecasts and temporally correlated trajectories, then clear the market with a two-stage linear program with recourse actions, as proposed in [START_REF] Morales | Pricing electricity in pools with wind producers[END_REF]. Here, the DA schedule considers the forecast uncertainty and anticipates possible network congestion. • Stochastic market clearing with PF (PF-Stoch): Derive predictive prescriptions by embedding the two-stage linear program [START_REF] Morales | Pricing electricity in pools with wind producers[END_REF] within the proposed tree algorithm. The key difference with the FO-Stoch is that scenarios are now weighted based on prescriptive importance. On a high level, the two-stage linear program with recourse actions is formulated as

min z, y(ξ) E ξ   t∈[T ] c(z t , y t (ξ); ξ|γ = γ)   s.t.
DA and RT Nodal Balance ∀ξ ∈ Ξ, t ∈ [T ],

DA and RT Flow Limits ∀ξ ∈ Ξ, t ∈ [T ],

Generator Technical Limits ∀ξ ∈ Ξ, t ∈ [T ], (9)
where z t are the first-stage (DA) decisions, and y t (ξ) the recourse (RT) decisions that depend on the realization of uncertainty ξ, for each period t. The detailed mathematical formulation is provided in Appendix A.

For simplicity, we further assume that reserve capacity procurement does not incur additional costs, therefore, all flexible operators with spare capacity are eligible for providing upward and downward regulation during RT operation (w.r.t to technical constraints). As the stochastic market clearing formulation co-optimizes energy and reserves, this mainly affects the deterministic formulation, in order to avoid clearing an additional reserve capacity market. Finally, we highlight that our solution explicitly accounts for recourse actions, and network and temporal constraints when generating predictive prescriptions, which are not considered in relevant works [START_REF] Zhang | Cost-oriented load forecasting[END_REF], [START_REF] Morales | Prescribing net demand for electricity market clearing[END_REF].

2) Experiment setup: The modified IEEE-24 bus system described in [START_REF] Ordoudis | An updated version of the ieee rts 24-bus system for electricity market and power system operation studies[END_REF] is used in this experiment. Regarding load, we employ data from the Global Energy Forecasting Competition 2014 (GEFCom2014) [START_REF] Hong | Probabilistic energy forecasting: Global energy forecasting competition 2014 and beyond[END_REF]. A total of 2.5 years of load data is selected, with the last year reserved for testing. Features γ comprise day-ahead temperature forecasts (averaged over all meteorological stations), historical lags and relevant calendar variables. We assume a peak load of 2700MW and scale the GEFCom data between the range of [0.20 × 2700, 2700]. The load is distributed to the respective nodes according to the percentages provided in [START_REF] Ordoudis | An updated version of the ieee rts 24-bus system for electricity market and power system operation studies[END_REF]. The cost of load shedding is set at 500$/M W h. To avoid infeasible solutions during RT operation, an additional slack is included for downward regulation. The cost is set significantly lower than load shedding, at 200$/M W h. Regarding model training, the PF-Stoch hyperparameters are set at {B = 20, K = 5, n min = 10}, with individual trees trained in parallel. Scenario reduction is discussed when applicable.

3) Results: Table II presents the aggregated out-of-sample costs as a function of sample size n, with PF-Stoch resulting in the lowest overall costs. Further examination shows that the PF-Stoch leads to the highest expected DA costs and the lowest expected redispatch costs, as illustrated in Fig. 4. Similarly to the previous experiment, the difference between the PF-Stoch and the FO solutions is more pronounced for smaller n, with all models converging for 1.5 year of data. For a sample size of 6 months, the learning component of the FO methods fails to generate adequate load forecasts. The stochastic market clearing design (both FO-Stoch and PF-Stoch) leads to lower aggregated costs, as a result of reduced redispatch costs during RT operation which is to be expected, as the stochastic formulation endogenously models the forecast uncertainty. The deterministic market clearing (FO-Det) consistently results in lower DA costs, also to be expected.

To gain further insight, Fig. 5 examines a single day of operation for the full training data. Fig. 5a presents the aggregated DA production, which matches the respective load forecast, and the actual aggregated demand. Note that both approaches incur similar forecast error for the selected day. In Fig. 5b, the DA dispatch schedule of the flexible generators is presented for 7:00, 15:00, and 23:00. In all cases, the PF-Stoch approach results in a larger number of generators being scheduled, with fewer of them reaching their maximum capacity. In turn, this leads to higher availability of regulation services, in both directions, during RT operation. Interestingly, this holds true even for 7:00 and 23:00, where the aggregated scheduled production for the PF-Stoch is significantly smaller than for the FO-Det and FO-Stoch approaches. This is attributed to the PF assigning scenario weights according to their prescriptive importance. In contrast, in the FO-Stoch solution all scenarios are equiprobable; therefore, a scenario that could possibly lead to large redispatch costs is considered solely based on probability of occurrence, which mitigates its impact on downstream costs. Overall, we conclude that the proposed modeling approach improves the utility of the stochastic market clearing design. Finally, we examine the impact of scenario reduction during learning. The results presented in Table III are obtained using a standard PC featuring Intel Core i7 CPU with a clock rate of 2.7 GHz and 16GB of RAM (no parallelization), for the full dataset, and hyperparameters {B = 20, K = 5, n min = 10}. We observe that a higher reduction in CPU time correlates with a larger performance decrease, with a maximum reduction of around 81% in CPU time for n max = 10. For n max = 50, a 43% reduction in CPU time is associated with a modest 2.87% decrease in performance. Note that increasing K, i.e., the number of splits evaluated per node, is expected to result in even larger relative reductions in CPU time. Overall, the embedded scenario reduction provides an attractive method to improve scalability with a modest decrease in performance.

IV. CONCLUSIONS

This work presented a data-driven method for value-oriented forecasting and decision-making in the presence of contextual information. The proposed method is based on prescriptive trees, trained under alternative loss function, that approximate the solution to conditional stochastic optimization problems. Color density ranks generators according to cost (higher density denotes higher cost). Fig. 5: DA schedule for a single day.

An internal scenario reduction methodology to facilitate training is also described. Experiments in two case studies examine the performance compared to deterministic and stochastic lookahead policies derived under the standard "forecast, then optimize" approach. The results show that the proposed method improves the out-of-sample prescriptive performance, and is less sensitive to the training sample size. Overall, the results validate the proposed approach as an effective alternative to the standard modeling approach. Future work could focus on adapting the proposed approach in an online setting and tackling sequential decision-making problems.

(RT) decisions are further augmented by index ξ to refer to a specific scenario.

1) Deterministic formulation: The deterministic formulation of the market clearing problem with network constraints (multi-period DC-OPF) is as follows

minimize Z t∈[T ] c p t (10a) subject to A G p t -A L d t = Bθ da t ∀t ∈ [T ], (10b) 0 ≤ p t ≤ p ∀t ∈ [T ], (10c)
-R d ≤ p t -p t-1 ≤ R u ∀t ∈ [T ], (10d)
-f ≤ b diag Aθ da t ≤ f ∀t ∈ [T ], (10e)
θ da,ref t = 0 ∀t ∈ [T ], (10f)
where Z = {p, θ da } the set of DA decisions and d t the expected nodal demand (forecasts). The objective (10a) minimizes the DA dispatch costs, subject to nodal balance constraints (10b), generator technical limits (10c), ramping constraints (10d), and line capacity limits (10e), with (10f) setting the reference node.

2) Redispatch problem: When uncertainty is realized during RT operation, the operator achieves the demand-supply equilibrium by solving the redispatch problem. Here, the optimal solutions (p

-f ≤ b diag Aθ rt ξ t ≤ f ∀t ∈ [T ], (11f)
θ rt,ref ξ t = 0 ∀t ∈ [T ], (11g)
where Ψ = {r u ξ t , r d ξ t , s u ξ t , s d ξ t , θ rt ξ t } defines the set RT decisions under scenario ξ . Specifically, r u , r d ∈ R |G| denote the generators upward/downward regulation, s u ∈ R |L| denotes a slack variable for upward regulation (i.e., load shedding), s d ∈ R |G| denotes a slack variable for downward regulation, and θ rt ∈ R |V| denotes the RT voltage angles. The RT redispatch problem minimizes the linear balancing costs (11a), subject to RT nodal balance constraints (11b), technical constraints for the corrective actions (11c)-(11f), and RT flows (11g), (11h). Slack s d is used to avoid infeasible solutions, by imposing an additional cost for downward regulation exceeding the maximum downward reserve capacity R -.

3) Stochastic Market Clearing: : stochastic formulation [START_REF] Morales | Pricing electricity in pools with wind producers[END_REF] explicitly models recourse actions, thus reserves are estimated endogenously. For ξ ∈ Ξ, where Ξ = {ξ i } n i=1 is a finite set of scenarios, the two-stage linear program is defined as

The predictive prescriptions generated with the proposed PF approach are based on the same formulation, with the objective (12a) being modified according to the learned weights ω n,i (x), as follows

Algorithm 2

 2 ReducedSAA Input: Data {ξ i } n0 i=1 , hyperparameter n max Output: SAA solution of subproblem with scenario reduction 1: Apply scenario reduction, generate new set of samples ξ ∈ Ξ , where |Ξ | = n max 2: Estimate cost v = min z∈Z 1 |Ξ | c(z; ξ i ) and prescription z = arg min z∈Z 1 |Ξ | c(z; ξ i ) 3: Scale cost to aggregate level v aggr = n 0 v 4: return (v aggr , z)

Fig. 1 :

 1 Fig. 1: Example of DA load forecasts: point forecasts, probabilistic forecasts as prediction intervals (PI), and scenarios.

Fig. 2 :Fig. 3 :

 23 Fig. 2: Plots of aggregated profit versus penalties for the various sample sizes. Marker size is analogous to the value of γ, . Points towards the top and left are better.

Fig. 4 :

 4 Fig. 4: Expected daily operating cost.

  DA dispatch schedule of flexible generators for specific hours.

  ) -(10f) ∀t ∈ [T ] (12b) constraints (11b) -(11h) ∀ξ ∈ Ξ, t ∈ [T ]

  min Z, Ψ ξ t∈[T ]   c p t + ω n,i (x) ξ∈Ξ c r u ξt , r d ξt , s u ξt , s d ξt , s u ξt   (13a) s.t. constraints (10b) -(10f) ∀t ∈ [T ],(13b) constraints(11b) -(11h) ∀ξ ∈ Ξ, t ∈ [T ].

TABLE I :

 I Problem Parameters

	Parameter	Value
	z soc	1
	z in	0.5
	z out	0.2
	z soc , z in , z out	0
	η in	0.8
	η out	0.9
	γ,	{0.01, 0.05, 0.1, 0.5, 1}

TABLE II :

 II Aggregated Operational Cost (10 6 $)

			Sample Size n FO-Det FO-Stoch PF-Stoch
			6 months	119.39	115.69	80.04
			1 year	83.49	81.90	73.84
			1.5 years	82.20	81.70	73.77
		600		FO-Det	FO-Stoch	PF-Stoch
		500			DA Costs DA Costs	RT Costs RT Costs
	$)
	Expected Cost (10 3	200 300 400
		100
		0	6 months	1 year Sample Size n	1.5 years

TABLE III :

 III Impact of Scenario Reduction on Average CPU

	Total Scheduled Production (MW)	1000 1200 1400 1600 1800	Actual FO-Det/Stoch PF-Stoch
		0	5	10	Hour	15	20
	(a) Aggregated DA dispatch schedule. Both approaches incur a Mean
	Absolute Error of 116.1MW for the selected day.
							Time and Performance
							PF	PF with Scenario Reduction, n max
							10	20	50
							CPU time (min)	18.05 ±1.42 3.35 ±0.38 5.74 ±0.62 10.26 ±1.2
							Cost increase (%) -	17.71%	7.46%	2.87%

  * t , θ da * t ) define problem parameters. For a specific realization ξ of uncertainty d, the problem is formulated as follows

	minimize Ψ	t∈[T ]	c r u ξ t , r d ξ t , s u ξ ,t , s d ξ t	(11a)
	subject to
	A G (r u ξ t -r d ξ t -s d ξ t )
	+ A (11c)
	0 ≤ r d ξ t ≤ min(p * t , R	-) ∀t ∈ [T ],	(11d)
	0 ≤ s d ξ t ≤ p * t -r d ξ t ∀t ∈ [T ],	(11e)
	0 ≤ s u ξ t ≤ d ξ t ∀t ∈ [T ],

L (s u ξ t -d ξ t + d t ) = B(θ rt ξ t -θ da * t ) ∀t ∈ [T ], (11b) 0 ≤ r u ξ t ≤ min(p -p * t , R + ) ∀t ∈ [T ],

APPENDIX

A. Mathematical formulation for the market clearing problem

Let G denote the set of generators, L the set of inelastic loads, Λ the set of transmission lines, and V the set of buses. We further define B ∈ R |V|×|V| to be the bus susceptance matrix, A ∈ R |Λ|×|V| to be the bus-branch incidence matrix, b diag ∈ R |Λ|×|Λ| to be a diagonal matrix with line susceptances, and A G ∈ R |V|×|G| and A L ∈ R |V|×|L| to be auxiliary matrices that map generators and loads to buses. The generation p ∈ R |G| and the voltage angles θ da ∈ R |V| define the first-stage (DA) decision variables, while the uncertain demand is denoted as d ∈ R |L| . All symbols are augmented by index t to refer to a specific period. The forecast horizon is T = 24, and [T ] is shorthand for {1, . . . , T }. Second-stage


Douglas Rossi Melo

Cesar Albenes Zeferino
email: zeferino@univali.br

Eduardo Augusto Bezerra
email: eduardo.bezerra@ufsc.br

Luigi Dilillo
email: dilillo@lirmm.fr

Design and Evaluation of Implementation Impact on a Fault-Tolerant Network-on-Chip Router

Keywords: SoC, NoC, Fault Tolerance

This work investigates synthesis alternatives to minimize error propagation in the controllers responsible for flow regulation, packet routing, and resource arbitration in a Network-on-Chip router. The controllers are based on Finite-State Machines to provide flexibility and favor low resource usage in programmable logic devices. The proposed router embeds hardening techniques by using triple modular redundancy on controllers and the Hamming code on buffers. Experimental results show that the packet routing controller has the most impact on the metrics evaluated and that the migration from a Moore to a Mealy controller implementation reduces the error propagation and offers a higher throughput than hardening the controllers. The main contribution of this work consists of assessing the impact of different implementations of a router in terms of error propagation.

I. INTRODUCTION

Reducing the component size and increasing the operating frequency of integrated circuits makes the Systemson-Chip (SoCs) more susceptible to faults. Thus, depending on its operational environment, an SoC requires faulttolerant components to minimize error propagation [START_REF] Sorin | Fault tolerant computer architecture[END_REF], including the communication infrastructure. As current multiand many-core systems use Networks-on-Chip (NoCs) as interconnection architecture, the NoC must be able to detect a fault and prevent the resulting error from causing an application failure. However, providing reliability in an NoC affects silicon costs, communication performance, and power consumption, as it is usually done through redundancy.

Studies about fault tolerance in NoCs mainly address both transient and permanent faults. For instance, [START_REF] Zimmer | A fault model notation and error-control scheme for switch-to-switch buses in a network-on-chip[END_REF]- [START_REF] Pereira | Mechanisms to provide fault tolerance to a network-on-chip[END_REF] examine Single Event Upset (SEU) in NoC components. The studies [START_REF] Zimmer | A fault model notation and error-control scheme for switch-to-switch buses in a network-on-chip[END_REF]- [START_REF] Park | Exploring fault-tolerant network-on-chip architectures[END_REF] address transient faults that are due to crosstalk. The authors of [START_REF] Cota | A high-fault-coverage approach for the test of data, control and handshake interconnects in mesh networks-on-chip[END_REF]- [START_REF] Chen | Towards maximum utilization of remained bandwidth in defected noc links[END_REF] investigate the problem of short and open circuit faults in the links of a Networkon-Chip. In works [START_REF] Loi | Characterization and implementation of fault-tolerant vertical links for 3-d networks-onchip[END_REF]- [START_REF] Eghbal | Analytical fault tolerance assessment and metrics for tsv-based 3d network-on-chip[END_REF], the authors discuss the yield of vertical links in 3-D NoCs. The work [START_REF] Pinheiro | Optimized fault-tolerant buffer design for network-on-chip applications[END_REF] proposes a fault-tolerant buffer design. The studies [START_REF] Cota | A high-fault-coverage approach for the test of data, control and handshake interconnects in mesh networks-on-chip[END_REF], [START_REF] Pasricha | A low overhead fault tolerant routing scheme for 3d networks-on-chip[END_REF], and [START_REF] Feng | Addressing transient and permanent faults in noc with efficient fault-tolerant deflection router[END_REF] examine intermittent faults and treat them as permanent faults. All of these works address the use of techniques to protect a given component and do not explore the hardening of different components using different techniques.

In this context, this work aims at evaluating the performance and resilience of an NoC router using combinations of flow regulation (or flow control), routing, and arbitration controllers, presenting the trade-off between the use of hardware resources and the susceptibility to error propagation. The proposed router implements hardening techniques through the application of the Triple Modular Redundancy (TMR) on controllers and the Hamming code on buffers. The first versions of router architecture were presented in [START_REF] Melo | Analyzing the error propagation in a parameterizable network-on-chip router[END_REF] and [START_REF]Maximizing the inner resilience of a network-on-chip through router controllers design[END_REF], where we explored the combination of different router controllers.

In this work, we protect the router's buffers and combine it with the controllers protection for evaluation in terms of silicon cost, communication performance, and reliability. The main contribution of this work is the assessment of the impact of different implementations on the inner resilience of the router itself in terms of error propagation.

The remainder of this paper is structured in the following sessions. Session II describes the architecture of the proposed fault-tolerant router, while Session III presents the methods to evaluate the router. Following, Session IV presents and discusses the experimental results, and Session V gives the final remarks.

II. ROUTER ARCHITECTURE

A. The Baseline Router

This work presents a router architecture for 2-D mesh networks. The router has five ports named Local, North, East, South, and West. The Local port is the terminal at which a processing core is attached, and the other ports connect the router to its neighbors. Internally, it uses a distributed architecture with input and output channels interconnected by a crossbar (Fig. 1). The input channels comprise controllers responsible for input flow regulation and packet routing, while the output channels include controllers, which perform the channel arbitration and output flow regulation. Each input channel also integrates a parameterizable input buffer. The flow regulation controller implements a 4-stage handshake protocol for receiving and sending packet flits (a flit is the smallest piece of data over which is performed the flow regulation). The routing controller runs the XY algorithm to request an output channel to forward an incoming packet. The arbitration controller consists of a Round-Robin arbiter that schedules the use of the output channel by the packets in the input channels. All these controllers are composed of Finite-State Machines (FSMs) and the router was designed to enable each controller to use a different type of FSM (Moore or Mealy). This approach enables us to better investigate architectural trade-off.

Input Channels

B. The Hardened Router

The proposed router architecture implements fault tolerance techniques. Since the memorization elements are the most susceptible to SEU faults, the internal controllers and the input buffers of the router were hardened.

The Triple Modular Redundancy (TMR) was the technique chosen to protect the FSM of each controller as it can mask an error transparently [START_REF] Sorin | Fault tolerant computer architecture[END_REF]. It consists of replicating the component in three units operating over the same input signals and using a major voter to compare each controller's output and elect the most common output value, as illustrated in Figure 2. This technique implies a high resource overhead if applied to complex structures. However, in the case of components with few interface signals and memory elements, as the controllers used in the proposed router, the TMR technique can provide a good cost/benefit. Hamming code [START_REF] Hamming | Error detecting and error correcting codes[END_REF] was used to protect the input buffers of the router. This technique relies on interleaving parity bits to detect up to two errors and correct a single error in a data word. Its implementation consists of an encoder, which generates the parity bits that compose the Error-Correcting Code (ECC) field, and a decoder that performs the error correction. We use Hamming solely to correct a single error. However, it is possible to forward a double error detection to an external structure, such as a network interface, for retransmission requests.

Controller

The encoder used for the Hamming code technique comprises a set of XOR gates, while the decoder checks parities through a Look-Up Table . In case of a mismatching output, the decoder fixes the wrong bit by inverting its value. As the ECC information must be attached to each packet flit, the number of memory elements of the input buffers increases. Figure 3 shows the application of the Hamming code.

Encoder Decoder Input Data

Output Data

III. MATERIALS AND METHODS

We described the baseline and the hardened routers using VHDL. Afterward, we synthesized them to FPGA using the Intel Quartus Prime design suite for costs and performance assessment. We selected the Intel Cyclone V 5CGTFD9E5F35C7 as the target device and disabled all the synthesis tool's optimizations options. Costs express the usage of resources-i.e., Look-Up Tables (LUTs) and Flip-Flops (FFs)-reported by the synthesis tool. Performance is given by the maximum operating frequency and the maximum throughput. The former is obtained by using TimeQuest Analyzer, and the latter is computed by dividing the number of bits transferred through all the communication channels by the simulated time.

Figure 4 shows the packet format used for router verification. It consists of a single bit to perform frame control, a single flit as the header, two payload Flow Control Units (flits), and a trailer. The header flit solely addresses the coordinates of the destination router. Both the header and the last payload flit (trailer) use '1' as the frame bit, while the regular payload flits use '0'.

For resilience assessment, we simulated a fault injection campaign on a single router using Mentor Graphics ModelSim simulator. For assessment, we defined a worst-case scenario with a workload that continuously injects 4-flit packets (i.e., injection rate equals 100%) to a fixed set of non-concurrent routes inside the router, thus enabling an assessment to be done at maximum loading, as shown in Fig. 5. Fault coverage is obtained by the number of propagated errors normalized in space, considering the amount of FFs, and along the time, adjusted by the throughput. For fault injection, we adopted the strategy proposed by the authors of [START_REF] Travessini | Processor core profiling for seu effect analysis[END_REF], which injected SEU faults into the registers of a processor using built-in commands of the ModelSim simulator, causing a bit-flip. Each iteration of the fault injection strategy includes the following steps:

1) Simulation without fault injection to obtain a golden run.

2) Listing of all registers in the circuit and picking of a random one to inject a fault.

3) Randomly determination of the time occurrence of the fault within the simulation time window, 4) Simulation pre-fault injection. 5) Fault injection: bit flip forced into the selected register. 6) Simulation for a predefined time interval. 7) Comparison between the outputs outcoming the fault injection and the golden runs.

In the experiments, we evaluated the baseline and the hardened router in configurations without any fault-tolerant mechanism (which we refer to as STD -standard) or applying TMR to the controllers or the Hamming code to the buffers (which we refer to as HAM). We also varied the types of FSM used to implement the flow regulation (input and output), routing, and arbitration controller, resulting in eight different architectural configurations for each router.

The configurations used 32-bit wide flits and buffers capable of storing up to 4 flits. Each configuration performed 1,000 simulations running for 100 µs (i.e., the simulated time). This approach was applied to obtain a more accurate measurement of the error propagation rate among all scenarios.

IV. RESULTS

A. Baseline

Table I presents the synthesis results for the baseline router. Note that all the controllers and buffers use the socalled standard (STD) implementation. As expected, we see that most of the Mealy-based controllers require more LUTs, due to the additional decoding of the output signals in states, and require less FFs because the FSMs encode few states. Moreover, increasing combinational logic implies a longer critical path and a lower operating frequency. Thus, the maximum operating frequency of the fully Moore approach is 70% higher than that of the configuration using only Mealy controllers. We highlight that the five buffers require the use of 895 LUTs, on average, and 695 FFs, being responsible for at least 50% of the combinational area and 90% of the sequential area, depending on the configuration used. This aspect has impact on the effectiveness of each fault-tolerant technique, as we show below.

B. Hardened Controllers

Table II presents the synthesis results for the hardened controllers using TMR (the buffers are not protected). The technique increased the number of logical elements needed, which reflected in a longer critical path and consequent performance degradation. In comparison with the baseline router, we see an average increase of 16.9% in combinational logic and 12.2% in sequential logic. The maximum operating frequency decreases 14.2%, on average. Fig. 6 presents the maximum throughput for the different configurations of the FSM type used in the baseline (STD) and hardened (TMR) implementations. Throughput is computed considering each configuration running for 10 Kcycles at its maximum operating frequency. As we can see, the higher the operating frequency, the highest the throughput. Moreover, we see that the use of TMR degrades performance in 13.7%, on average, and the highest levels of throughput are obtained by using the Mealy machine in the flow regulation controller and the Moore FSM in the arbitration controller. Fig. 7 presents the error propagation rate for the different configurations. We see that the use of TMR to protect the controllers has low effectiveness in reducing error propagation. Moreover, in configurations with an error propagation rate higher than 15%, the use of TMR propagates even more errors than the unprotected version. This result is because the faults affect the router randomly, and most of the faults are injected into the buffers' registers, which are not protected by the TMR technique. After, we show how reliability can be improved by hardening these buffers.

C. Hardened Buffers

Table III presents the synthesis results for a router configuration in which only the buffers are hardened (i.e., the controllers are not protected). Buffer protection resulted in a similar overhead in combinational (average of 533 LUTs) and sequential (120 FFs) logic resources among all combinations. The maximum operating frequency decreases when the Mealy machine is used for routing.

Figure 8 presents the throughput of the standard and hardened buffers configurations. We can see that the use of Hamming code on input buffers imposed a performance degradation of approximately 40% in comparison with the baseline router configuration. This high reduction observed in the maximum operating frequency is due to the increase in the critical path, given by the Hamming decoder at the buffers' output and the combinational logic of the Mealy-based routing controllers. Regarding the error rate, presented in Fig. 9, we see that the use of the Hamming code technique effectively reduces error propagation. As discussed earlier, these results confirm that the buffers are the components most affected by the injected faults.

D. Hardened Controllers and Buffers

As a last experiment, we analyzed the indicators for the fault-tolerant router implementing the TMR technique on the controllers and Hamming code on the buffers. Table IV presents the costs in terms of the resources occupied in all the configurations. As expected, these implementations have the highest occupancy of logic resources. As the router area increases, the longer wires increase the critical path and degrades performance. The maximum operating frequency and the throughput (Fig. 10) of this router are even lower than those of the configurations that implement only one fault-tolerant technique. Despite the high silicon overhead inferred by the hardening techniques, the router's resilience significantly increases in this configuration when comparing with the baseline router. As we can see in Fig. 11, the error propagation of a full-hardened router is less than 2%, on average.

E. Discussion

The experimental results enable designers and system integrators to understand some aspects related to the impact of the architectures of the controllers and the fault-tolerant techniques in the indicators of cost, performance, and reliability of the router, as follows:

• The adoption of the Mealy machine in the controllers increases the router's performance when applied to the flow controllers and resilience when applied to the routing controllers, because of the reduced number of states in the FSM compared to Moore implementations. • Although the TMR technique is admittedly expensive, the impact of its use in controllers is not significant since most of the router's cost is due to buffers. On the other hand, the hardening of the controllers alone has low effectiveness in reducing the error propagation precisely because most of the faults are on the buffers. • The Hamming code technique is highly effective in reducing the propagation of errors. However, this technique significantly degrades the router's performance, mainly when the routing controller is based on the Mealy machine because this combination results in a longer critical path than those of the other configurations. Despite presenting a cost and performance degradation, the use of Hamming code is justified in applications reliability. Configurations that use Mealy in the flow regulation and Moore in the routing controller are those that present the highest throughput. When configured with a Mealy FSM for arbitration, it offers the lowest error propagation.

Although being designed to handle transient faults, the developed router may tolerate permanent single-bit stuck faults since their effect is covered by the embedded fault tolerance techniques. Furthermore, the router provides some protection against multi-bit faults, as long as they occur on different controllers and different buffer flits. However, a physical test campaign is required to obtain realistic fault coverage for given environment such as atmosphere or space.

V. CONCLUSION

This paper presented the design and assessment of a fault-tolerant router for an NoC. The experimental results enable a system designer to select the configurations of the NoC that best fit the application's requirements and costs. For example, an on-board satellite computer requires high reliability in communication, while signal processing applications prioritize throughput.

As future work, we plan to extend the verification campaign to analyze the error propagation rate when injecting multiple faults. We also intend to conduct the physical assessment in a particle accelerator to deepen the reliability analysis. The source code of the proposed architecture is available at [START_REF] Xarc | [END_REF].

Fig. 1 .

 1 Fig. 1. Proposed router architecture.

Fig. 2 .

 2 Fig. 2. Triple Modular Redundancy on controllers.

Fig. 3 .

 3 Fig. 3. Hamming code on buffers.

Fig. 4 .

 4 Fig. 4. Packet format for verification.

Fig. 5 .

 5 Fig. 5. Traffic scenario.

Fig. 6 .

 6 Fig.6. Throughput of standard and hardened controllers (Gbit/s).

Fig. 7 .

 7 Fig. 7. Error rate of standard and hardened controllers.

  Fig. 8. Throughput of standard and hardened buffers (Gbit/s).

Fig. 9 .

 9 Fig. 9. Error rate of standard and hardened buffers.

  Fig.10. Throughput of standard and hardened controllers and buffers (Gbit/s).

Fig. 11 .

 11 Fig. 11. Error rate of standard and hardened controllers and buffers

  Throughput of standard and hardened buffers (Gbit/s).

		22.0%	23.1%
					15.9%
	15.0%
	10.0%			9.0%
	5.0%	5.2%	3.6%	4.8%	2.6%	5.1% 1.8%	6.2%	5.4% 3.0%
	0.0%

  Throughput of standard and hardened controllers and buffers (Gbit/s).

	22.0%	23.1%
			15.9%
		10.4%
		9.0%
				5.1%	6.2%
	2.37%	2.09%	1.44%	1.18%	2.04%

ACKNOWLEDGMENTS

This work was financed in part by CNPq -the Brazilian National Council for Scientific and Technological Development -Processes 315287/2018-7 and 436982/2018-8.