Ying Tang
email: ying.tang@univ-lorraine.fr

Alessio Franci
email: afranci@ciencias.unam.mx

Romain Postoyan
email: romain.postoyan@univ-lorraine.fr

Romain Postoyan Qualitative

Qualitative parameter estimation for a class of relaxation oscillators

Keywords: parameter estimation, relaxation oscillator, singular perturbation, neuroscience, nonlinear system, Lyapunov method

  L'archive ouverte pluridisciplinaire

Qualitative parameter estimation for a class of relaxation oscillators

Ying Tang * , * * Alessio Franci * * * Romain Postoyan * , * *

INTRODUCTION

Online parameter estimation of dynamical models of neuronal activity might lead to new perspectives in neurosciences. In epilepsy for instance, having access in real-time to the gains governing the excitation/inhibition balance within populations of neurons might provide important information about the on-going electrophysiological activity, and thus help developing new strategies to detect or even predict seizures. This task is challenging because of the redundancy and the large variability of biophysical parameters across populations of neurons and neuronal circuits exhibiting similar activity patterns. Disparate combinations of biophysical parameters are indeed known to lead to the same activity pattern at the cellular level [START_REF] Goldman | Global structure, robustness, and modulation of neuronal models[END_REF]. The same degenerated parametrization property propagates at the neuronal circuit level [START_REF] Marder | Variability, compensation, and modulation in neurons and circuits[END_REF]. In addition, biophysical parameters slowly vary, inducing sharp transitions between qualitatively different activity modes (spiking or bursting, healthy or epileptic, etc.) at the crossing of critical parameter sets. In this context, quantitative modeling and fitting of experimental data generically constitute an ill-posed problem and a new estimation approach is therefore needed.

Online parameter estimation for dynamical systems has been widely studied, in particular for linear systems, see

This work was supported by the ANR under the grant SEPICOT (ANR 12 JS03 004 01), by the French "Région Grand Est" through a fellowship grant 2016-2017, and by DGAPA-UNAM under the grant PAPIIT IA105816. e.g., [START_REF] Ioannou | Robust adaptive control[END_REF]; [START_REF] Ljung | System identification[END_REF]. Online estimation technique for nonlinear systems is still a developing field and the available results are more scarce, see e.g., [START_REF] Bastin | Stable adaptive observers for nonlinear time-varying systems[END_REF]; [START_REF] Farza | Adaptive observers for nonlinearly parameterized class of nonlinear systems[END_REF]; [START_REF] Adetola | Finite-time parameter estimation in adaptive control of nonlinear systems[END_REF]; [START_REF] Besançon | Nonlinear observers and applications[END_REF]; [START_REF] Mauroy | Linear identification of nonlinear systems: A lifting technique based on the Koopman operator[END_REF]; [START_REF] Moles | Parameter estimation in biochemical pathways: A comparison of global optimization methods[END_REF]. Recently, results have appeared in the context of neurosciences. For instance, in [START_REF] Schiff | Kalman filter control of a model of spatiotemporal cortical dynamics[END_REF], the authors have applied unscented nonlinear Kalman filters to estimate the states and the parameters of a spatiotemporal model of the cortex. Other stochastic filtering methods such as Monte Carlo schemes and Bayesian approach have also been considered for computational neuroscience applications [START_REF] Brockwell | Recursive bayesian decoding of motor cortical signals by particle filtering[END_REF]; [START_REF] David | Dynamic causal modeling of evoked responses in EEG and MEG[END_REF]. While these methods can provide very interesting results, these also have a high computational cost in general, may be difficult to tune, and there is usually no convergence proof. On the other hand, deterministic schemes motivated by neuroscience applications have recently been provided in [START_REF] Tyukin | Adaptive observers and parametric identification for systems in non-canonical adaptive observer form[END_REF]; [START_REF] Chong | Parameter and state estimation of nonlinear systems using a multi-observer under the supervisory framework[END_REF], in which the estimation errors are ensured to converge to the origin under persistency of excitation conditions. All the estimation approaches reviewed in this introduction rely on some type of quantitative fitting of the measured data, which questions their applicability in the neuronal context.

In the present work, we cast the problem in a deterministic setting and we propose a different approach. Indeed, our objective is no longer to asymptotically estimate the unknown parameters, but first to identify analytically those which rule the type of behavior of the model, and second to estimate the distance between these parameters and critical values at which a change of activity occurs. We illustrate these ideas on a class of nonlinear systems with a single sigmoidal nonlinearity and two sharply separated timescales. In order to capture salient properties of neuronal dynamics, we suppose that the exact functional form of the nonlinearity is unknown, as well as the timescale separation. We firstly show that, independently of the particular expression of the sigmoid nonlinearity, a single parameter rules the transition between global asymptotic stability and almost global convergence to a relaxation limit cycle. This part of the paper is inspired by the work in [START_REF] Drion | Dynamic input conductances shape neuronal spiking[END_REF]. The parameter redundancy has algebraically been tackled by extracting the few ruling parameters governing neuronal dynamics and their critical values. Those parameters define qualitative models that provide a constructive geometric framework to analyze modulation and robustness of neuronal activity in a principled way Franci et al. (2014). Based on the assumption that the input and the states of the system are known, we subsequently design an online qualitative estimator that estimates the distance of the ruling parameter to the unknown critical value at which the rest/oscillation transition happens. The designed estimator therefore provides online information about the system activity type and whether a change of activity is prone to appear without any quantitative fitting of the measured data. The estimator design is guided by singularity theory [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF] and its convergence and robustness properties are ensured by the joint use of Lyapunov analysis and geometric singular perturbation theory [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF]. The idea of qualitative estimation is also related to the problem of steering a system toward an a priori unknown bifurcation point [START_REF] Moreau | Balancing at the border of instability[END_REF].

The proofs of the results are omitted for space reasons, and these can be found in [START_REF] Tang | Qualitative parameter estimation for a class of relaxation oscillators[END_REF].

Notation. The usual Euclidean distance is denoted by | • |. For a function h : R 0 → R, the associated infinity norm is denoted by h [0,∞) = sup s∈[0,∞) |h(s)|, when it is well-defined. We use sgn(x) to denote the sign function from R to {-1, 0, 1} with sgn(0) = 0. For any function f : R → R, we denote the range of f as f (R) = {z : z = f (x) for some x ∈ R}. Let A, B be two nonempty subsets in R n , their Hausdorff distance is noted by

d H (A, B) = max{sup a∈A inf b∈B |a -b|, sup b∈B inf a∈A |a - b|}. A continuous function γ : R 0 → R 0 is of class K ∞ if it
is zero at zero, strictly increasing and unbounded.

PROBLEM STATEMENT

A class of two-time-scale nonlinear systems

We consider the nonlinear system ẋf = -x f + S(βx f + u -x s ), ẋs = ε(x f -x s ), (1a) (1b) where x f , x s ∈ R are the state variables, β ∈ R is an unknown parameter, u ∈ R is the input, which is assumed to be a known constant, and 0 < ε 1 is a small parameter, which is unknown. The mapping S : R → R is an unknown sigmoid function, which is assumed to satisfy the following properties.

Assumption 1. The following properties hold. a) S is smooth. b) S(0) = 0. c) S (a) > 0 for all a ∈ R (monotonicity), and argmax a∈R S (a) = 0 (sector-valued). d) sgn (S (a)) = -sgn(a) for all a ∈ R.

2

Standard sigmoid functions such as a → c 1 tanh(c 2 a), a → c 1 arctan(c 2 a), a → c1 1+e -c 2 a -c1 2 with suitable c 1 , c 2 ∈ R, verify Assumption 1.
Due to the small parameter ε, system (1) evolves according to two time scales. We follow the standard approach of singular perturbation theory [START_REF] Kokotović | Singular perturbation methods in control: analysis and design[END_REF]; [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF]) to decompose system (1) into two subsystems, which represent the fast dynamics and the slow dynamics, respectively, called the layer and the reduced subsystems. By setting ε = 0 in (1b), we obtain the layer dynamics behavior 2b), the slow variable x s is treated as a constant parameter.

ẋf = -x f + S(βx f + u -x s ), ẋs = 0. (2a) (2b) In (
To account for slow variations of x s , we rescale time as

τ = εt, hence d dt = ε d dτ . Then system (1) becomes εx f := -x f + S(βx f + u -x s ), x s = x f -x s , (3a) (3b)
where stands for d dτ . Setting ε = 0 in (3a), we obtain the reduced dynamics

0 = -x f + S(βx f + u -x s ), x s = x f -x s . (4a) (4b)
The reduced dynamics evolves in the slow time τ and is an algebro-differential equation. In particular it defines a one-dimensional vector field over the critical manifold

S 0 u,β := {(x s , x f ) ∈ R 2 : -x f + S(βx f + u -x s ) = 0}.
(5) The critical manifold S 0 u,β depends on u and β. However, since we assume u and β constant, for simplicity we omit the index u and β in the rest of the paper. The critical manifold will therefore simply be denoted by S 0 . Let

β c := 1 S (0) , (6)
which is well-defined according to item c) of Assumption 1. This parameter will play a key role in the sequel.

Stability analysis

We first state stability properties for system (1) in the next proposition.

Proposition 1. Consider system (1), the following holds.

(1) For any β -β c < 0 and any constant input u, there exists ε > 0 such that, for all ε ∈ (0, ε], system (1) has a uniformly globally exponentially stable fixed point. Moreover, all trajectories converge in an O(ε)-time to an O(ε)-neighborhood of the critical manifold S 0 defined in (5).

(2) a) For all 0 < β -β c < 1, there exists ū > 0, such that, for any constant input u in (-ū, ū), there exists ε > 0 such that, for all ε ∈ (0, ε], system (1) possesses an almost globally attractive and locally exponentially stable periodic orbit P ε . In particular, for any initial condition, except that in the unique unstable fixed point of (1), the solution converges to P ε . Moreover, P ε is strongly attractive, i.e. any solution in a neighborhood of the periodic orbit P ε converges to it exponentially fast with decay rate e -K/ε , for some constant K > 0. b) Let T ε be the period of P ε and p : R → P ε , with p(t) = p(t + T ε ) for all t ∈ R, be the associated periodic solution of (1). Let P 0 be the singular limit of P ε . The Hausdorff distance between P ε and P 0 verifies d H (P ε , P 0 ) = O(ε 2/3 ). Moreover, for all t 0 ∈ R, there exists δT ε ⊂ [t 0 , t 0 + T ε ) and

|δT ε | = O(ε)T ε , such that |p(t) -S 0 | < O(ε) for all t ∈ [t 0 , t 0 + T ε ) \ δT ε . 2
Proposition 1 indicates that the solution to system (1) converges either to a stable fixed point or to a stable limit cycle. The transition between the two behaviors depends only on the single parameter β -β c and not on the exact expression of the sigmoid S.

Remark 1. The proof of Proposition 1 relies on differential geometry arguments ( [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF], Theorem 1,2,3), (Krupa and Szmolyan, 2001b,a, Theorem 2.1)). We mention that the approach in [START_REF] Kokotović | Singular perturbation methods in control: analysis and design[END_REF]; [START_REF] Khalil | Nonlinear systems (3rd Edition)[END_REF] for singularly perturbed nonlinear systems is not applicable in the context of this work because, when β > β c , the critical manifold of system (1) possesses singularities, which implies that (4a) no longer has isolated roots as often assumed. 2

Objective

Our objective is to detect online in which type of activity system (1) is, that is, whether solutions converge to a fixed point or a limit cycle (oscillation) and whether we are close to a change of activity. According to Proposition 1, the two types of behaviors depend only on the value of β.

In the simple, academic setting of this paper, we might approximate the sigmoid function S by Taylor expansion up to some order and quantitatively estimate the expansion coefficients together with ε and β. However, due to the generic impossibility of quantitatively modeling and fitting experimental neuronal data, this approach would not overcome, in a real experimental setting, the limitation of quantitative estimation reviewed in the introduction. Instead, we aim at qualitatively estimating the distance of the ruling parameter β from critical value without using any type of quantitative fitting. Indeed, this distance is all we need to assess online whether the system exhibits oscillation or has a globally exponentially stable fixed point and whether it is near to a change of activity.

Proposition 1 proves that almost all trajectories of system (1) converge to an O(ε)-neighborhood of the critical manifold S 0 in both the fixed point and the limit cycle cases, at least for most time in the latter case. Based on singularity theory [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF], Franci and Sepulchre (2014) and in view of Assumption 1, the qualitative shape of the critical manifold S 0 is the same as that of the set

{(x f , x s ) : -x 3 f + (β -β c )x f + u -x s = 0}, ( 7
) whose definition is independent of S. Qualitatively, the shape of this set, as well as of S 0 , is fully determine by the sign of β -β c (see Fig. 1). When β < β c , it is the graph of a monotone decreasing function. For β = β c , it has a point of infinite slope, corresponding to the hysteresis singularity, and, when β > β c , it is S-shaped. We exploit this information in the next section to design and analyze the estimator.

QUALITATIVE ESTIMATOR

Estimator design

In the following, we construct the estimate of β -β c based on the sole information provided by ( 7) about the qualitative shape of the critical manifold. We first make the following assumption. Assumption 2. Both the fast x f and the slow x s variables are measured in system (1) and the control input u is known and constant.

2

We propose the following nonlinear parameter estimator β = -kx f (-

x 3 f + βx f + u -x s ) := f ( β, x f , u -x s ), (8)
where k > 0 is a design parameter. We note that f is independent of S and ε, its expression is only related to the set (7). The influence of the design parameter k and the robustness of the estimator will be briefly discussed in Section 4.

Steady-state properties of the qualitative estimator

Steady-state of β satisfies -x 3 f + βx f + u -x s = 0. Hence, we implicitly define the function β * (x f , x s , u) such that

-x 3 f + β * (x f , x s , u)x f + u -x s = 0, (9)
for u, β ∈ R and x f , x s on S 0 , i.e. x f , x s are related by the critical manifold equation

-x f + S(βx f + u -x s ) = 0. ( 10
)
Recalling that the sigmoid function S is invertible on S(R) according to item c) of Assumption 1, we can explicitly invert (10) as follows

x s = -S -1 (x f ) + βx f + u. (11)
Replacing this expression for x s in (9), we obtain the explicit expression of β * (x f , x s , u)

β * (x f , x s , u) = x 3 f -S -1 (x f ) + βx f x f . ( 12
)
We see that β * depends only on x f and β. We therefore write β * (x f , β) in the following.

The next lemma provides important properties of β * (x f , β).

In particular, it ensures that β * (x f , β) is well-defined when x f = 0. Furthermore it formally characterises in which sense β * provides qualitative information about β -β c . Lemma 1. The following holds. 1) β * is smooth on S(R) × R.

2) For all β ∈ R,

∂ β * (x f ,β) ∂β
> 0 when β < β c and

x f ∈ S(R), while ∂ β * (x f ,β) ∂β > 0 when 0 < β -β c < 1 and x f is in a neighborhood of the origin. 3) β * (x f , β) = β -β c + O(x 2 f ) for any x f ∈ S(R) and β ∈ R.
2

Lemma 1 is informative about the qualitative dependence of β * on β -β c . According to item 2) of Lemma 1, β * is strictly increasing in β, and hence in β -β c , in the resting activity (i.e. β < β c ). If a change of activity is prone to appear (i.e. β is close to β c ), then necessarily x f will lie in a neighborhood of the origin (for u small), which relates to the point of infinite slope, corresponding to the hysteresis singularity of the critical manifold. In this case, β * will change its sign from negative to positive according to item 3) of Lemma 1, and then it will grow according to the second part of item 2) of Lemma 1. The same holds in the other direction, meaning that if β * goes from positive to negative values. Hence, β * allows to detect a change of activity in system (1). If β converges to β * as time tends to infinity, then we can distinguish the different behaviors of the system based on β. This is the purpose of the next subsection.

Stability analysis

Firstly, we discuss stability property of ( 8) with (x f , x s , u) as input in Lemma 2. Then we study the different behaviors of the three-dimensional overall system, i.e. system (1) and estimator (8) in Proposition 2. Finally, we show the link between β and β * in Theorem 1.

The next lemma states an incremental input-to-state stability property (in the semiglobal sense) of system ( 8) with (x f , x s , u) as input. In the sequel we use the triplet (x f , x s , u) to denote a solution to system (1) with arbitrary initial condition and input u. Lemma 2. Let ∆ > 0 and S ∆ be the set of solutions (x f , x s , u) to system (1) such that: (i) max x f [0,+∞) , x s [0,+∞) , u [0,+∞) ∆; (ii) there exist t * (∆) 0 and a(∆) > 0, such that for any t t * (∆),

t t * (∆) x 2 f (τ ) dτ a(∆) (t -t * (∆)). ( 13
)
Then, for (x f 1 , x s1 , u 1 ), (x f 2 , x s2 , u 2 ) ∈ S ∆ , for any M > 0 and | β1 (0)|, | β2 (0)| M , there exist strictly positive constants ∆,M , ω ∆,M and K ∞ function γ ∆,M , such that the corresponding solutions β1 , β2 to system (8) satisfy for any t 0

| β1 (t) -β2 (t)| ∆,M e -ω ∆,M t | β1 (0) -β2 (0)| +γ ∆,M x f 1 -x f 2 [0,t) + x s1 -x s2 [0,t) + u 1 -u 2 [0,t) . ( 14
)
2

Similar to semiglobal ISS with respect to inputs as defined in [START_REF] Angeli | Further equivalences and semiglobal versions of integral input to state stability[END_REF], property ( 14) indicates that the estimator ( 8) is semiglobal incremental ISS with respect to initial conditions and inputs. Condition ( 13) is similar to a persistency of excitation condition. In the present paper, we assume it and we verify it online in simulations.

Conditions on the initial states of system (1) and its input to ensure it a priori will be investigated in future work.

To finalize the analysis, we consider the overall system, that is (1) and ( 8), which is three-dimensional and given by

     β = -kx f (-x 3 f + βx f + u -x s ), ẋf = -x f + S(βx f + u -x s ), ẋs = ε(x f -x s ). (15a) (15b) (15c)
Note that system (15) is a cascade of system (15b)-( 15c) with (15a). We state the convergence property of system (15) in the next proposition. Proposition 2. Consider system (15).

(1) For any β -β c < 0, there exists ε > 0 such that for any ε ∈ (0, ε], x f (0), x s (0) ∈ R and constant input u such that the corresponding solution (x f , x s , u) to (15b)-(15c) is in the set S ∆ for some ∆ > 0, and for any β(0) ∈ R, the corresponding solution to system (15) asymptotically converges to a fixed point. (2) For any 0 < β -β c < 1, there exists ε > 0 such that for any ε ∈ (0, ε], x f (0), x s (0) ∈ R and constant input u such that the corresponding solution (x f , x s , u) to (15b)-(15c) is in the set S ∆ for some ∆ > 0, and for any β(0) ∈ R, the corresponding solution to system (15) asymptotically converge to a periodic limit cycle. 2

Proposition 2 ensures that the estimate β converges either to a fixed point or to a limit cycle depending on the sign of β -β c . Now we are ready to state the link between the estimate β and β * defined by ( 12) in the following theorem. Theorem 1. Consider system (15).

(1) For any β -β c < 0, let u and ε be defined as in item (1) of Proposition 2, and let x * f be the x f -component of the fixed point of (15). Then β(t) converges to β * (x * f , β) as time goes to infinity. (2) For any 0 < β -β c < 1, let u and ε be defined as in item (2) of Proposition 2, and let x lc f (t + θ) be the T ε -periodic asymptotic orbit of x f (t) along the limit cycle to which converges x f . Then, for all t 0 0, there exists

δT ε ⊂ [t 0 , t 0 + T ε ) with |δT ε | = O(ε)T ε , such that β(t) -β * (x lc f (t + θ), β) = O(ε) for all t ∈ [t 0 , t 0 + T ε ) \ δT ε . 2
When β < β c , according to item 1) of Proposition 2, the estimate converges to a fixed point as time increases, and item 1) of Theorem 1 specifies that its value is β * . When 0 < β -β c < 1, the estimate converges to a periodic limit cycle in view of item 2) of Proposition 2, and item 2) of Theorem 1 indicates that the estimate value is in an O(ε)neighborhood of β * for all the time, after a sufficiently long time, except during jumps and each is of length

|δT ε | = O(ε)T ε .
We see that β * is not a constant in this case but a state-dependent signal, which is unusual in parameter estimation. In practice, we can average the obtained value of β over a sliding window.

NUMERICAL ILLUSTRATIONS

We analyze the results provided by the proposed qualitative estimator in simulations in this section. We first consider the sigmoid function S : x → tanh(x), for which β c = 1. We choose ε = 0.001, u = -0.01. Figure 2 illustrates the time evolutions of the states x f , x s of system (1) for different values of β -β c . We find that x f converges to a constant value as time increases, as well as x s , when β -β c < 0, which is consistent with item 1) of Proposition 1. Both x f and x s converge to an oscillatory behavior when 0 < β -β c < 1, which is in agreement with item 2) of Proposition 1.

We have designed the estimator as in ( 8) with k = 5.

To check that the estimator is able to detect a change of activity of the system, we use a step signal for β.

In particular, β = 0.9 on [0, 500] and β = 1.1 on [500,1000], which leads to a change of sign of β -β c at t = 500. We observe in Fig. 3 that, when β -β c < 0, the estimate converges to a constant value as time increases, corresponding to the resting activity. When 0 < ββ c < 1, the estimate tends to a periodic function, which is strictly positive. Hence, it indicates the oscillation activity. Moreover, the estimate value evolves from negative to positive when the change of sign of β -β c occurs. This implies the change of the system's behavior. It could be interpreted as, for example in epilepsy a change from rest to seizure. We have tested different values of k in (8). The simulations indicate that the speed of convergence of β increases with k. However, this also leads to bigger spikes at "jumps", which may provide negative estimate value during a very short interval.

To test the robustness of our approach, we next consider the following sigmoid functions:

S 1 : x → tanh(x), S 2 : x → x √ 1+x 2 , S 3 : x → 1 1+e -x -1 2 .
The corresponding critical values are β c1 = β c2 = 1 and β c3 = 4. We emphasize that even though the nonlinearities are different, the estimator remains the same as defined in (8). Figure 4 shows the relationship between β and β -β c for different input u and perturbation parameter ε. Note that, when β -β c < 0, the value of β is chosen as the constant value, to which the estimate converges. When 0 < β -β c < 1, the estimate tends to a periodic behavior. Then β is selected as the average value of this periodic function (over a period). We observe that for β-β c < 0, the estimate value β is negative and it increases as β -β c increases. When 0 < β -β c < 1, β is positive. Moreover, when β -β c = 0, the estimate is around the origin. Hence the estimate β qualitatively infers the distance of β -β c to the origin. We remark that the estimator works robustly to variations of the sigmoid shape.

To further evaluate the robustness of the scheme, we have added small additive measurement noises given by d out (t) = 0.005 sin(50t) for x f , x s and d in = 0.008 for input u. It is found from Fig. 5 that the estimator still provides good results in this case.

CONCLUSION

We have, to the best of our knowledge, introduced for the first time the concept of qualitative estimation, which allows to detect online the behavioral state of a multipletimescale nonlinear system, independently of large uncertainties on the system nonlinearities and without using any type of quantitative fitting. This is achieved by first extracting the system ruling parameter(s) and by subsequently designing an estimator of the distance of the ruling parameter(s) to some critical value(s) at which the behavioral transition happens. Timescale separation is instrumental in the estimator design, because it allows to reduce the qualitative estimation problem to the estimation of the qualitative shape of the system critical manifold. As a first illustration, we have focused on a class of twodimensional nonlinear systems with two time-scales and a single nonlinearity either exhibiting resting or relaxation oscillation behaviors. Future extensions will focus on generalizing the theory to higher-dimensional systems, different unknown nonlinearities, and multiple possible qualitative behaviors. We will also investigate the case where only the fast variable is measured, and not the full state as done in the paper.

Fig. 1 :

 1 Fig. 1: Geometrical form of the critical manifold

βFig. 2 :Fig. 3 :

 23 Fig. 2: States x f and xs of system (1) for different values of β.

Fig. 4 :Fig. 5 :

 45 Fig.4: Robustness with respect to different sigmoid functions. Vertical axis is β and horizontal axis is β -βc. Different S functions:S 1 : x → tanh(x) (red solid line), S 2 : x → x √ 1+x 2 (blue dashed line), S 3 : x → 1 1+e -x -12 (pink - * line).


Kristian Pagh Nielsen

Frank Vignola

Lourdes Ramírez
email: lourdes.ramirez@ciemat.esd

Philippe Blanc
email: philippe.blanc@mines-paristech.fre

Richard Meyer
email: richard.meyer@suntrace.def

Manuel Blanco
email: m.blanco@cyi.ac.cy

Excerpts from the Report: "BeyondTMY -Meteorological Data Sets for CSP/STE Performance Simulations"

In order to facilitate comprehensive economic modeling of CSP/STE power plants realistic long-term meteorological datasets with temporal resolution down to 1 minute is a main premise. Currently available standard datasets do not fulfil this premise. The datasets also need to combine the high quality of well-maintained ground-based irradiance measurements and the global coverage of satellite-derived data. Even with the best available data it is necessary to account for the uncertainty in this and the sampling uncertainty from finite time-series to enable the optimal statistical characterization. It is a general challenge that satellite-derived data lack the required temporal resolution, and also often does not cover periods with major volcanic eruptions. Here we see prospects in synthetically generated realistic datasets, although research and development work is required on how to optimally produce and quality assure these.

INTRODUCTION

To build a solar thermal electric (STE) power plant an investment of the order of hundreds of millions euro is needed. An investment of that size requires bank support. The banks need to have an idea of the plant production for several years into the future. A main variable here is the future variability of the meteorological conditions, where direct normal irradiance (DNI) is the primary variable of importance for concentrating solar power (CSP). This enables estimations of the project profitability (the financial gain to be expected) and the pay-back period (the time required to recover an investment) [1]. The data for financing purposes should represent as well as possible the local conditions to be expected at the plant site at least over the tenure of the loan, which typically is in the order of 15 years. Some investors interested in the 'golden end' of a power project might also consider up to 20-25 years. Thus, in case of CSP/STE projects: 1) The long term average of solar resources should well represent the average; 2) The interannual variability of the DNI should be represented; 3) The short-term variability of DNI, which also plays a major role on yields of CSP/STE-plants [START_REF] Chhatbar | The Influence of Meteorological Parameters on the Energy Yield of Solar Thermal Power Plants[END_REF], should include realistic data sequences and interdependencies between the meteorological variables.

Yearly datasets of hourly solar data [START_REF] Benseman | [END_REF], meteorological reference years [START_REF] Andersen | Referenceåret -Vejrdata for VVS beregninger, (The Reference Year -Weather data for HVAC-calculations)[END_REF][START_REF] Lund | The "Reference Year", a set of climatic data for Environmental Engineering[END_REF][START_REF] Lund | Design Reference Years[END_REF], typical meteorological years (TMYs) [START_REF] Hall | Generation of Typical Meteorological Years for 26 SOLMET Stations[END_REF][START_REF] Marion | User's Manual for TMY2s[END_REF][START_REF] Wilcox | User's Manual for TMY3 Data Sets[END_REF], all consist of 12 typical months based on a long-term (at least 10 years) meteorological data set.

These datasets have been a main data source for simulations for decades. They are, however, not optimized with respect to the DNI solar resource and they lack untypical months.

The yearly datasets also do not represent the interannual variability. To accommodate this lack, they can be supplemented with yearly probabilities of exceedance (PoE) values that estimate the solar resource in bad years. The PoE is the probability that a given level of yearly power generation will be exceeded. Thus P90 is the yearly production level of a system that is expected to be exceeded in 90% of all years and P99 is the yearly production expected to be exceeded in 99% of all years [START_REF] Vignola | Statistical Analysis of a Solar-Radiation Dataset: Characteristics and Requirements for P50, P90, and P99 Evaluation[END_REF][START_REF] Cebecauer | [END_REF]. Consensus is needed on how to derive these PoE values. Here we show and discuss examples of the challenges in deriving them.

The temporal resolution in the yearly datasets is mostly 1 hour. One exception is the Design Reference Year (DRY) developed in the IEA SHC Task 9, in which DNI data with 5 minute resolution is recommended [START_REF] Lund | Design Reference Years[END_REF]. Since the short-term variability of DNI plays a major role for CSP/STE plant energy yields, optimal meteorological datasets should be developed that include 1 minute resolution DNI data. The difference between hourly and 1 minute resolution DNI data on two days with variable cloud cover is illustrated in Fig. 1. In recent years, the growing number of STE projects has pushed researchers to look for solutions to specific needs of this technology. To generate TMY or reference yearly data sets several meteorological variables are weighted and compared on a monthly basis to their long term average distributions. 12 specific months are then chosen based on the fit of the weighted variables (e.g. [START_REF] Hall | Generation of Typical Meteorological Years for 26 SOLMET Stations[END_REF][START_REF] Marion | User's Manual for TMY2s[END_REF]). Since DNI is the single most important variable for CSP/STE yield, new yearly datasets with 12 specific months chosen only from the fit of DNI have been developed [START_REF] Hoyer-Klick | Characteristic meteorological years from ground and satellite data[END_REF][START_REF] Habte | Temporal and Spatial Comparison of Gridded TMY, TDY, and TGY Data Sets[END_REF]. In the ENDORSE TMY generation service [START_REF] Espinar | Report on the Product S4: 'TMY for production[END_REF], the maximum wind gust threshold can also be taken into account. The mirrors in a CSP/STE plant need to be placed in safe mode during strong wind gusts.

Another important development in recent years is the production of gridded data sets (e.g. [START_REF] Hoyer-Klick | Characteristic meteorological years from ground and satellite data[END_REF][START_REF] Habte | Temporal and Spatial Comparison of Gridded TMY, TDY, and TGY Data Sets[END_REF]). In these satellite-derived irradiances are combined with ground-based measurements to obtain global coverage that includes locations for which long-term ground-based datasets are not available. The satellite-derived irradiances need to be adapted to quality controlled ground-based measurements during at least one year and preferably two years. This process is referred to as site-adaptation [START_REF] Meyer | Towards Standardization of CSP Yield Assessments[END_REF][START_REF] Ramírez | Towards the standardization of procedures for solar radiation data series generation[END_REF] and has recently been reviewed by Polo et al. [START_REF] Polo | [END_REF].

The upcoming IEC/TC117 standard for STE plants will include a technical specification (TS 62862-1-2) on yearly data set generation for STE simulations. This specification is based on several years of research, but does not include currently discussed issues such as statistical distributions, uncertainty and variability. This is because industrial standards are based on consensus and include only agreed minimum requirements. In Fig. 2 the timeline of making such an industrial standard is illustrated. This illustrates why supplementary guidelines are needed that also include current state-of-the-art methods. The relation of the new guidelines to the current standards will be presented. For financing CSP/STE plants yearly datasets generated according to TS 62862-1-2 can be used as a base case. For detailed feasibility studies and profitability assessments additional characterizations of long-term variability and uncertainty are needed. The BeyondTMY project is aimed at making ambitious recommendations in this regard form the research and development forefront.

BeyondTMY is a SolarPACES project that was started in 2015. In the project, methods for creating meteorological data sets for CSP/STE performance simulations have been reviewed, and an updated roadmap for research, development and deployment is being developed. A comprehensive report is being made in 2016. The full report will be available at the SolarPACES 2016 conference. This paper includes excerpts and our main findings from this report. We briefly describe the main data sources, the uncertainty of these, sources of DNI variability that are challenging to account for, statistical characterization of long-term datasets and our outlook for the opportunities and requirements for synthetically generated datasets.

DATA SOURCES

Ground-Based Data

The most accurate DNI measurements are made with absolute cavity radiometers (ACR). The ACR DNI measurements are used as reference values because electrical current can be measured to a much higher degree of accuracy than measurements of the thermal energy that other instruments measure. The internationally recognized standard for DNI is the World Radiometric Reference (WRR) developed and maintained by the World Radiation Center [START_REF] Finsterle | WMO International Pyrheliometer Comparison, IPC-XI: final report[END_REF]. Periodically, other ACRs are calibrated against this standard and they achieve accuracies, at the 95% level of confidence, in the range of ±0.36% to 0.40% [START_REF] Reda | NREL Pyrheliometer Comparisons[END_REF]. These calibrated ACRs and then used to calibrate other DNI measuring instruments. All certifiable calibrations of DNI measuring instruments can trace their calibrations to the WRR.

ACR's are very expensive and are not intended for field operation. There are four ways to obtain DNI values from ground-based measurements in the field: 1) Field Pyrheliometers, 2) Rotating Shadowband Irradiometers, 3) Pyranometers with a shadow masks, 4) Calculations using Global Horizontal Irradiance (GHI) and Diffuse Horizontal Irradiance (DfHI). In most cases there are many varieties of each type of DNI instrument.

Most field pyrheliometers consist of thermopiles-based detectors at the end of a collimation tube that has a window covering the aperture. The aperture and the collimation tube provide an opening with a full angle of view between 5.0 and 5.7, where newer models of pyrheliometers have 5.0 full angle field of view. Note that concentrating solar technologies may not utilize this exact 5.0 field of view but may utilize a narrower or broader field of view.

Other ground-based data used in meteorological datasets include: GHI measured with pyranometers, temperature and relative humidity measured in 2 meters height. Wind and wind direction measured in 10 meters height. The measurement heights are standards of the World Meteorological Organization (WMO) [START_REF] Wmo | Guide to Meteorological Instruments and Methods of Observation[END_REF].

Satellite-Derived Data

There are a variety of ways to obtain estimates of surface solar irradiance values from satellites observations and other meteorological measurements. The basis for any modeling of solar irradiance data is the ability to accurately estimate the irradiance during cloudless periods. There are three basic ways to estimate surface irradiance from satellite images. The first is to use radiative transfer models and cloud cover information that is part of the reanalysis methodology used for climate studies. These are called physical models and this method has been used for the NASA/SSE irradiance database [START_REF] Sengupta | Best Practices Handbook for the Collection and Use of Solar Resource Data for Solar Energy Applications[END_REF]. The empirical approach is based on statistical and correlation methods to model the irradiance from satellite images using ground-based measurements to determine the parameters and validate the model. Currently the best models use a combination of both approaches to take advantage of the physical information being gathered from satellites [START_REF] Polo | [END_REF].

The 2015 IEA SHC Task 36 & 46 best practices report [START_REF] Sengupta | Best Practices Handbook for the Collection and Use of Solar Resource Data for Solar Energy Applications[END_REF] provides a comprehensive overview of both satellitederived and ground-based irradiance data sources.

IRRADIANCE UNCERTAINTY

All data, whether measured or modelled, is uncertain. The uncertainty is important to quantify and account for in both simulations and statistical analyses. Uncertainty should not be confused with the actual variability of the data.

To standardize the discussion of uncertainties and to set standards that provide guidelines for determining uncertainty and characterizing uncertainties, the Guide to Expressing Uncertainties in Measurements (GUM) was created [START_REF]ISO/IEC Guide 98-3: Uncertainty of measurement -Part 3: Guide to the expression of uncertainty in measurement (GUM:1995)[END_REF]. This document explains in detail the GUM terminology and explains how to perform uncertainty analysis using the GUM procedures. The GUM model starts by defining the "measurand", the quality that is being measured. The resulting measurement is an approximation or estimate of the measurand and a full description of the measurement includes the uncertainty of the measurement. In addition, other environment quantities that affect the measurement should be included. For example, the WRR calibrations are measurements of DNI made when DNI values >700 W/m2 under clear sky, stable weather conditions.

Traditionally, errors are viewed as having two components, a random and a systematic component. Random errors arrive from unpredictable or stochastic temporal and spatial variations of quantities that influence the measurements. Random errors can usually be reduced by increasing the number of measurements. The experimental standard deviation is a measure of the uncertainty of the mean resulting from random effects. Systematic errors arise from a recognized effect that influences the measurements. Error has the connotation of mistake or fault while uncertainty has the connotation of ambiguity. Error and uncertainty have often been used interchangeably, but under the GUM methodology uncertainty more precisely describes the situation. The systematic effect can be quantified and if it is significant in size, an adjustment factor, also called correction factor in the literature, can be devised and applied to compensate for the effect. The uncertainty in the adjustment factor is a measure of incomplete knowledge of the value required for the adjustment. Therefore one has uncertainties in the measurements, possible adjustments made to account for some systematic effects plus the uncertainties associated with the modeled adjustments. The terms error and uncertainty should be used precisely and care should be taken to distinguish them.

The following is a sampling of possible sources of uncertainty in measurements from the GUM document.

IRRADIANCE VARIABILITY

To be able to correctly describe the variability of DNI it is important to understand the root causes of this. Both natural and anthropogenic causes affect this variability. Several of these causes of variability are only well represented in long-term meteorological dataset. A good example is the recent 2014-2016 El-Niño in the Pacific. An El-Niño event of similar strength last occurred in 1997-1998. The variability of this is quantified with the Multivariate El-Niño Southern Oscillation (ENSO) index (MEI) [START_REF] Noaa | Earth System Research Laboratory: Multivariate ENSO Index (MEI)[END_REF] or the ENSO 3.4 index that is based on Pacific sea surface temperatures and is available from the National Oceanic and Atmospheric Administration (NOAA) [START_REF] Noaa | Physical Oceanography Division (PHOD) of AOML: El-Niño Southern Oscillation ENSO 3.4[END_REF]. The ENSO is a chaotic variation that affects the cloud conditions in and around the Pacific Ocean, and in lower latitude regions around the world. For variations in the North Atlantic region the North Atlantic Oscillation (NAO) index can be used to quantify large scale variability in the tracks of extratropical cyclones over time scales of weeks and months. The cyclone tracks affect cloud occurrence and thereby the regional solar resource [START_REF] Pozo-Vázquez | [END_REF]27].

The variability of global scale volcanic eruptions poses an even greater challenge to represent due to the low frequency of these. Data affected by large volcanic eruptions should as the Pinatubo 1991 eruption are omitted from TMYs and other yearly datasets. Likewise, P90 and similar probability of exceedance values often are estimated with a disclaimer that major volcanic eruptions are not considered. With the latest paleoclimatic data [28], assessments of the major volcanic eruption variability can be made. Sulphate aerosol deposits in Antarctic and Greenlandic ice cores can be related to the volcanic sulphate aerosols that are ejected into the stratosphere. These aerosols in the stratosphere have an e-folding decay time of about 1 year [29] and thus affect the solar resource significantly. Judging from the analysis of 2500 years of paleoclimatic data [28] the likelihood of a volcanic eruption ejecting more than 0.1 aerosol optical depth (AOD) sulphate aerosols into the stratosphere within a 10-year period is 43%.

Correspondingly, the likelihoods are 22%, 7% and 1.6% for AOD larger than 0.2, 0.5 and 1.0, respectively. To put this into context, the Pinatubo 1991 eruption is estimated to have caused a maximum stratospheric AOD of almost 0.2. These estimates rely on studies of the most recent major volcanic eruptions. Since some of the past ice core sulphate deposits are biased positively by proximal high latitude volcanos, the AOD estimates are likely to be overestimated. Work is ongoing with better isotope analysis of volcanic ice deposits that promises to enable recognition of individual volcano sources and the possibility of distinguishing tropospheric and sulphate [START_REF] Toohey | European Geosciences Union General Assembly[END_REF]. Thus, improved volcano risk assessments are likely to become available in the coming years.

Since the DNI resource is strongly affected by aerosols, scenarios of future trends of anthropogenic aerosol emissions are worth considering.

STATISTICAL CHARACTERIZATION OF LONG-TERM DATASETS

The long-term inter-annual DNI variability can be characterized with yearly probability of exceedance values. For example, the P90 value corresponds to the yearly DNI value that will be exceeded at least 90% of the years, during the lifetime of the project. It corresponds to complementary of the percentile: the P90 value of the yearly DNI corresponds to its 10th percentile. More generally Pxx can be used, where xx can be for instance 75% or 95%. The Pxx values depend on the statistical characterization of the yearly variability and whether the yearly data are independent and identically distributed (IID). Different Pxx values will be obtained depending on this characterization.

For illustration purpose, we have selected a very long-term (36 years; 1978-2013) and high quality in-situ measurements of DNI from the pyranometric station in Eugene, Oregon. This station belongs to the radiometric network UO SRML (University of Oregon Solar Radiation Monitoring Laboratory). In Fig. 3 results are shown for the PoE values P50, P75, P90 and P99 when three different statistical distributions are assumed: The Normal (or Gaussian), the Weibull, and the Gumbel distributions. Here it can be seen that the P99 estimates are more affected by which distribution is assumed than the lower PoE values.

The corresponding P90 of the yearly direct normal irradiation for the empirical data distribution is 1238 kWh/m 2 , which is within the 95% confidence intervals of all P90 values for the three statistical distributions but clearly larger than these. The short-coming of deriving PoE values from empirical data distributions is that these are strongly limited by the data available as noticed by Fernández Peruchena et al. [START_REF] Fernández Peruchena | [END_REF]. The estimated statistical distributions are also limited by the available data, and by their representativeness of the actual variability. Monte Carlo simulations can be used to quantify the sampling uncertainty for given statistical distributions. In Fig. 4 results are shown based on 3000 Monte Carlo runs consisting of three sets of 1000 random years generated from the Normal, Weibull and Gumbell distributions, respectively. with six tested methods: NORMFIT is the normal distribution, WBLFIT is the Weibull distribution, GBLFIT is theGumbell distribution, ECDF is the empirical distribution, KDE is the kernel density estimation and CLT is the central limit theorem.

This illustrates how the uncertainty of the estimated P90 values depends on the length of the time series available. Nevertheless, it is important to note that this uncertainty is only "theoretical". Firstly, the measurement uncertainty in the yearly DNI values should be accounted for, as this also affects the P90 estimations [32]. Secondly, this Monte-Carlo methodology is based on the very important assumption that the yearly DNI is an independent identically distributed process. In the example of the yearly DNI dataset from Eugene, Oregon, this assumption is clearly questionable, and we can suspect that it is not an exception as opposed to the rule. Using statistical tools based on assumptions that are not clearly verified for their applications could be risky. There is a need of advanced statistical methodology to be able to detect and to handle non-identically distributed and/or non-independent process.

As a perspective, one can consider Monte-Carlo Markov Chain and fractional Gaussian noise random generators to be able to account for observed long-range correlation in time series [33,34,[START_REF] Handschy | Reducing Wind Power Variability by Quantifying Geographic Diversity Length Scales[END_REF]. These techniques of "realistic" or plausible time series generation meant to reproduce observed inter-annual correlation, trends, long-range memory (Hurst exponent) may be efficient to generate large numbers of representative synthetic very long term historical DNI dataset to be used to assess the PoE values and their related uncertainties. We believe that stochastic economic analysis of CSP/STE plant projects makes it possible to incorporate both the uncertainty and variability associated with the input values into the analysis to the model. Modelling of the variability of the electricity generation of a CSP/STE plant associated with the variability of the solar resource and the rest of meteorological variables on the site is critical. Given this and good estimates of other uncertainties, stochastic economic analysis of CSP/STE plant projects is made possible. In Fig. 5 a flowchart of the financial modeling is shown that includes uncertainties and variability in all input data. These mostly contribute to the annual electricity generation via the variability in the solar resource data. Realistic data in the sense that they reproduce the natural variability both on the 1 minute time scale and on the 100 year time scales -to include effects such as major volcanic eruptions -are needed. Since high quality measurements at these extremes of the time spectrum are not available, realistic synthetically generated data should be a main focus for future investigations.

FIGURE 1 :

 1 FIGURE 1: Comparison of two days of DNI data on an hourly and a 1 minute scale. These measurements are from the Baseline Surface Radiation Network (BSRN) station in Carpentras, France. They have been measured with a Kipp & Zonen CH1 Pyrheliometer at a 1 second sampling rate averaged to 1 minute and 1 hour data [12].

FIGURE 2 :

 2 FIGURE 2: Activity timeline showing the progress from research scientist methods to industrial standards.

a)

  Incomplete definition of the measurand; b) Imperfect reaIization of the definition of the measurand; c) Nonrepresentative sampling -the sample measured may not represent the defined measurand; d) Inadequate knowledge of the effects of environmental conditions on the measurement or imperfect; e) Measurement of environmental conditions; f) Personal bias in reading analogue instruments; g) Finite instrument resolution or discrimination threshold; h) Inexact values of measurement standards and reference materials; i) Inexact values of constants and other parameters obtained from external sources and used in the datareduction algorithm; j) Approximations and assumptions incorporated in the measurement method and procedure; k) Variations in repeated observations of the measurand under apparently identical conditions.

FIGURE 3 :

 3 FIGURE 3: Estimations of P99, P90, P75 and P50 of the yearly DNI and the corresponding 95% confidence interval, with three fitted probability density function (PDFs), specifically the Normal, Weibull, and Gumbel PDFs.

FIGURE 4 :

 4 FIGURE 4: Median errors (P50) and confidence intervals (CI) for different levels (75%, 90% and 95%) of the estimation of P90with six tested methods: NORMFIT is the normal distribution, WBLFIT is the Weibull distribution, GBLFIT is theGumbell distribution, ECDF is the empirical distribution, KDE is the kernel density estimation and CLT is the central limit theorem.

FIGURE 5 :

 5 FIGURE 5: Stochastic economic feasibility analysis flowchart.

The GUM methodology breaks the uncertainties into two types: A and B. Both types are quantified by variances and/or standard deviations. The Type A evaluation is calculated from series of repeated observations For the Type B evaluations, the variance is evaluated using available knowledge, for example the characteristics of the measuring device. In the Tables 1-3 overviews of the sources of uncertainty for pyrheliometers, rotating shadowband radiometers and pyranometers with shadow masks are given, respectively.

ACKNOWLEDGEMENTS

We acknowledge support from SolarPACES for travel expenses during the BeyondTMY project. Additional funding from the Energy Development and Demonstration Program (EUDP) of the Danish Energy Agency, The Energy Trust of Oregon, The National Renewable Energy Laboratory, and Bonneville Power Administration is also acknowledged.


Hideki Mori
email: mori@citl.titech.ac.jp

Twenty-One Things to Do with Scratch in the Twenty-First Century Classroom

Introduction

Globally, programming education for children is receiving widespread attention as a part of twenty-first century skills. Learning programming is not just about learning to code, but it also helps children learn how to learn. However, unlike traditional classes where students mostly learn from teachers' teaching, learning through programing (in a broad sense, "making") needs to be approached differently. To make programming education more successful in elementary school classrooms, students must be self-learners, and moreover, teachers need to encourage students to be self-learners by acting as facilitators, mentors, or peers (and occasionally as teachers). Classrooms have many physical limitations, and school curricula, especially elementary school curricula, have even more limitations, such as time and subjects. How can we introduce programming and making to twenty-first century classrooms? Upon its publication, "Twenty things to do with a computer" (Papert & Solomon 1971) impacted both teachers and researchers. In this paper, I will share twenty-one things that I have been doing with Scratch and the tiny programmable computer "Cricket," invented by MIT Media Lab, in Japanese elementary school classrooms to maximize students' learning through making, in over 700 hours of classes for 15 years.

Things to Do with Scratch or Other Programming Tools in the Classroom

Making is an extremely playful activity for students. While it is most important to have fun while engaging in these activities, it is also necessary to have structures to ensure that students do not spend their time in class just by having fun. Below are the items that I consider when I design and run classes at elementary schools.

Become Professionals

To learn through making, it is necessary to become active learners rather than passive learners. Encouraging students to become professionals or experts, such as inventors, designers, and researchers, makes it easier for them to find personal meanings in activities. When students have their own meaningful purposes, it becomes easy to find their way back when they get lost while tinkering with and exploring activities.

Make Students Small Teachers and Teachers Big Students

Encouraging students to be small teachers and actively helping other students who are facing difficulty with making is a good way to run classes. In making activities, various problems arise, as many as the number of students. Solving problems by themselves is a learning for students, but when they are faced with problems that cannot be solved by themselves, adequate advice from experts, including teachers and other students who understand the problem well, is needed at an appropriate time. As only a limited number of students can be supported by a teacher, students teaching other students is an effective means to run classes.

Write Out a Goal in the Class

Creative Learning Spiral [START_REF] Resnick | All I Really Need to Know (About Creative Thinking) I Learned (By Studying How Children Learn) in Kindergarten[END_REF]) is an effective way to learn from making. To realize the Creative Learning Spiral process in a time-constrained class, it is necessary to understand the activities of "Imagine" -"Create" -"Play" -"Share" -"Reflect." By writing out the purpose of these activities in the class, students can direct the activities within the class and create a chance for "Imagine."

Work Collaboratively with Other Students

Working on activities through collaboration with multiple people is one of the fundamental approaches for constructionist learning to create opportunities to learn through making. In addition, by encouraging individual works and collaborative works, teachers can prevent students from relying on groups without doing anything on their own, and students can engage in effective interaction between individual learning and group learning.

Reflect on Activities in the Class

Reflecting on the activities of each class is done in the latter part of the class and is divided into two parts. In the first part, students reflect on the activities of the class and summarize them in their notebook so that they can share their findings with the class. In the second part, students listen to other students' findings and summarize their reflections again, including what they want to try the next time.

Share Findings in the Class

To share what students discovered in the class, a short presentation time at the end of the class works well for students. Also, to create informal sharing opportunities during the class, it is important for students to be able to walk around the classroom when they find something interesting. Students are always active in sharing ideas informally, and it is often found that students make similar project to the ones sitting nearby.

Be Good Imitators

Students may not like imitating or being imitated. Teachers need to make them understand that they can imitate good ideas and incorporate them into their work actively. While doing so, teachers also need to encourage students to come up with new ideas and respect the student who created the original idea.

Present Own Projects

It is important to make time to present students' projects to other students. If the making class is held over several number of classes, it is good to have enough time to prepare for presentations in addition to short presentations in each class. By summarizing thoughts in a presentation and receiving opinions from other students, students are likely to get new ideas.

Evaluate Oneself by Oneself

Unlike other classes, it is difficult to uniformly determine learning goals in classes where students work on their own making projects. Therefore, it is important to set each student's own goal by himself or herself, and evaluate the activities on his or her own in reflection parts.

Make Failures Meaningful

Debugging in programming is the most powerful tool for students' learning (Papers 1980). Most students do not like facing failures. Teachers need to tell students that making mistakes in the class is good, and it is most important for them to find their own answers from their mistakes.

Give Chances to Re-Try

To create opportunities for students to learn from trial and error, teachers need to allocate more time to students for programming and making things on their own. In some cases, it is impossible to try new ideas unless the work is fundamentally changed. By giving students opportunities to work on more than two times, it becomes easier to newly incorporate the experience of the previous project and apply ideas gained from other students.

Ask Students to Help in Preparing the Classes

Students are likely to participate actively in classes when they themselves participate in the preparation for classes instead of simply working in prepared environments. It is good to ask students to gather materials for their making projects and bring them to classes. It is also beneficial to assign students to prepare the class with teachers (although it sometimes takes more time for teachers to prepare the classes with students.)

Ask Students to Clean up the Classroom

It is hard for teachers to take part and clean up everything when using a lot of materials in tasks such as robot making. When students do the cleaning by themselves, it becomes an opportunity for them to reflect on the making process as they disassemble and clean up the projects they made.

Sharpen the Progress of the Class

When students concentrate on making, it often becomes difficult to shift to the next activities, such as listening to teachers. It is good to ask students to listen to teachers when teachers talk. However, it is still and always difficult for them to move from their making activity to the next activity. Instead of forcing students to follow orders like a computer program does, teachers need to understand how students actually concentrate on making.

Become Non-Teaching Teachers

Teachers cannot teach students everything in the class of making, especially when new technologies are introduced in the class. By giving up teaching everything and trying to explore unknown things together with students, teachers themselves are released from the fear of dealing with new educational materials.

Look at the Form of Students' Elementary Learning

In subject study, many students have learned related knowledge at the cram school in Japan, so it is difficult to see the elementary picture of students. Making activities also provide an opportunity to see the figure facing students' activities to check how they face learning itself. (However, recently, there are more places outside the school where students experience programming and robot making, so it is necessary to pay attention.)

Leave Students for a While

Students' interests are diverse. Even if it seems for teachers to be far from students' goals, guiding students may be a necessary process for them. In order not to teach them too much, teachers need to leave students for a certain period of time. However, technical problems need to be supported immediately.

Think with the Body

When students engage in programming, they only think with their brain. It is important for students to learn to use not only their brain but also their body when they are learning [START_REF] Papert | Mindstorms: Children, Computers, and Powerful Ideas[END_REF]. It is a good introduction for students that they become like computers and move their body according to commands, such as "raise your hand" when teachers ask (program) them to do so.

Think about Computers and Programming Related to Students' Familiar Things

Computers are embedded in a lot of things around students, from appliances to toys and video games, and these things are controlled by computer programs. Thinking about how their familiar things work and making those things by themselves also deepen their understanding of computers and programming.

Tell Students Programming Is for Learning How to Learn

Many students, parents, and teachers still think that programming education is necessary for mastering programming (coding) skills. Teachers need to tell them that programming is useful for learning how to learn, and making things as well as programming are the same ways of learning how to learn. (I had an experience where a first-grade girl told me that programming is unnecessary for her as she wants to be a florist in the future.)

Think about How to Make Use of the Experiences in Other Classes

Both students and teachers can use their making experiences in classes for other subjects. Thinking about how students and teachers can make use of them concretely, it is possible to introduce new learning into ordinary classes as well, and it will also provide a chance to add one more thing to the nextgeneration (twenty-second century) classroom.

Acknowledgments

I would like to thank all the students and teachers who work with me in classes. I would especially like to thank Manabu Sugisawa as most of the ideas in this paper emanated from the classes that I conducted with him. This work was supported by JSPS KAKENHI Grant Number 17K01110.


Gaili Gao

Dongxing Zhang

Jun Liu

Design of a New Soil-Tuber Separation Device on Potato Harvesters

Keywords: Potato harvester, Soil-tuber separation, Separation device, dynamics analysis

A new soil-tuber separation device on potato harvesters has been developed. It mainly comprises a linear motor, a rectangular screen, a frame, four supporting springs and four dampers. It is the first time to propose that the soil-tuber separation device is driven by the linear motor. Forces acting on the rectangular screen are found. A differential equation for describing the forces is set up. An equation for evaluating exciting force of the linear motor has been derived from the solution to the differential equation. An equation of the factors affecting on the device resonance has also been obtained. Additionally, an equation for evaluating the deformed magnitude of the damping springs in the dampers has been derived. And finally a motion simulation experiment for the rectangular screen on a potato harvester has been done with the Matlab/Simulink module. The experiment expressed that the method is feasible.

Introduction

Many methods were used to separate soil from potatoes. However, field separating operation caused most damage to potatoes. [START_REF] Curl | Potato harvester[END_REF] invented a tuber harvester having a primary bed and a secondary bed. The secondary bed is encompassed by a devining chain [START_REF] Curl | Potato harvester[END_REF]. [START_REF] Hyde | Potato harvester performance with automatic chain-load control[END_REF] suggested this type of tuber harvester could be equipped with a system for controlling the depth of material on the primary chain by automatically adjusting chain speed. The system provides uniform flow of material from the primary to the secondary chain over a wide variety of field conditions and ground speeds [START_REF] Hyde | Potato harvester performance with automatic chain-load control[END_REF]. Another similar potato harvester was described by [START_REF] Welp | Potato harvester[END_REF]. It comprises a screening conveyor belt and an elevator disposed after the screening conveyor belt. The elevator is formed from a double conveyor having an internally rotating bar-type conveyor belt with narrow openings and an externally rotating bar-type conveyor belt with wide openings [START_REF] Welp | Potato harvester[END_REF]. In order to improve soil elimination, shaking of the conveyor chain (bed) in a vertical or horizontal direction was applied to the potato harvesters [START_REF] Woodruff | A preliminary analysis of a high-frequency soil riddling device for use on potato harvesters[END_REF] [START_REF] Woodruff | A preliminary analysis of a high-frequency soil riddling device for use on potato harvesters[END_REF]. In addition, a pressure belt mechanism operating on the surface of the primary conveyor was proposed (Misener et al., 1989) [START_REF] Misener | Improved design for soil separation and vine removal mechanisms on a potato harvester[END_REF]. However, due to friction between chain (bed) or soil and potatoes, the potatoes were greatly damaged. A simple mechanical separator was developed by [START_REF] Feller | Development of a clod separator for potato packing houses[END_REF] . The machine utilizes the difference in the restitution coefficient between potatoes and soil clods.Both potatoes and clods are dropped from a feed conveyor onto a revolving steel roller where the resulting trajectories cause the potatoes to end up at a more remote location than the clods. The separator can not be used to separate stones from the potatoes [START_REF] Feller | Development of a clod separator for potato packing houses[END_REF]. Gan-Mor et al.(1986) used a appropriate elastic plate for reducing the restitution coefficient of stones without changing that of the potatoes or clods so that the stones were separated [START_REF] Gan-Mor | Absorbing stone impact to enable separation from potatoes[END_REF]. However, the impact between the steel roller or plate and potatoes caused damage to potatoes. [START_REF] Bouillon | Potato harvesting apparatus[END_REF] proposed another potato harvesting apparatus. Its separator has three primary cylindrical cages, followed by three secondary cylindrical cages and then three rubber-coated cylindrical cages. The separator can remove most of the dirt and break up the dirt clods [START_REF] Bouillon | Potato harvesting apparatus[END_REF]. However, much of mechanical injury was caused by the conveyor equipment. In addition, great power required to drive the separator was consumed. Misener and McLeod (1989) also tested a mechanism which consists of a combination of disc-shaped clod rollers and cylindrical brushes. The mechanism essentially eliminated the loose soil and, on the average, 60% of the stones and 47% of the clods [START_REF] Misener | A resource efficient approach to potato-stone-clod separation[END_REF].

The objective of this research is to develop a new soil-tuber separation device on potato harvesters. The device can reduce potato damage, reduce the power required to separate soil from potatoes, and improve the separation efficiency.

Specific work is divided into three parts:

(1) Model of soil-tuber separation device.

(2) Dynamics analysis for soil-tuber separation device.

(3) A motion simulation experiment for a rectangular screen. A model of a soil-tuber separation device on a potato harvester is shown in Fig. 1.

Model of Soil-Tuber Separation Device

The device mainly comprises a linear motor, a rectangular screen, a frame, four supporting springs ( two of the supporting springs aren't shown in Fig. 1), and four dampers. The rectangular screen consists of equidistantly spaced bars at the center and a connecting plate on each side of the bars. The frame has four inside slideways and two bottom supporting plates. The motionless part of the linear motor is held onto the frame by connecting components ( the connecting components aren't shown in Fig. 1). The motion part of the linear motor is attached to the connecting plates of the rectangular screen by another connecting components( the connecting components aren't shown in Fig. 1). The rectangular screen is connected to the supporting plates of the frame by the four supporting springs. The frame is equipped with the four dampers each of which has a spring, a guide, and a spring cover( the spring cover isn't shown in Fig. 1 ). The angle between the datum plane of the device and the level is β . When the device works the rectangular screen driven by the motion part of the linear motor vibrates up and down in inclined direction along the four slideways of the frame and the four guides of the dampers .

3 Dynamics Analysis for Soil-Tuber Separation Device

Evaluation of Exciting Force

The forces exerted on the rectangular screen in vibrating direction are (Fig. 2):

exciting force of the linear motor acting on the top of the rectangular screen, ; supporting force of the supporting springs acting on the bottom of the rectangular screen, ; damping due to surface friction from the motion of the rectangular screen, ; and component of gravitational force of the rectangular screen that includes the tuber and soil on it, Using vibrating direction as x axis and putting the stagnation point O at equilibrium position (Fig. 2), the supporting force in the x direction may be described as follows:

F N F ′ β cos mg
) (

s x iK N λ + = . ( 1
)
Where is number of the supporting springs, i K is stiffness of the supporting springs ( ),

1 m N s λ is the deformed magnitude of the supporting springs at equilibrium position( ), m x is the distance of the rectangular screen from the equilibrium position ( ).

m

The exciting force is defined by:

t F F ω sin 0 = . ( 2
)
Where is amplitude of the exciting force ( ), 0 F N ω is angular frequency of the exciting force ( ), is the exciting time ( ).

1 -

s rad t s

The damping may be expressed as follows [START_REF] Zhang | Mechanical Vibration[END_REF]:

x F e & γ = ′ . ( 3
)
Where e γ is damping coefficient (

).

m s N

Based on Newton's law of motion the differential equation for describing the forces acting on the rectangular screen is:

x m F N mg F & & = ′ - - + β cos . ( 4
)
Where is mass ( ), m kg g is acceleration of gravity (

), 2 s m β is the angle between the datum plane of the device and the level ( .).

deg At equilibrium position

s K i mg λ β = cos . ( 5
)
Substituting equation ( 1),( 2),(3),and(5) into equation ( 4) and rearranging terms gives:

t F x K i x x m e ω γ sin 0 = + + & & & If we define m F q and m K i p m n e 0 2 , 2 = = = γ then equation (6) becomes: t q x p x n x ω sin 2 2 = + + & & & . ( 7
)
A particular solution to equation ( 7) is:

) ( sin ) ( ψ ω - = t B t x . ( 8
)
where

2 2 2 2 2 4 ) ( ω ω n p q B + - = . ( 9
)
and

2 2 2 ω ω ψ - = p n tg . ( 10
)
Rearranging equation ( 9) yields:

2 2 2 2 2 1 1 ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ⋅ = p p n p p q B ω ω . ( 11
)
If we define p n and p K i

F p q B = = = = ζ ω λ , 0 2 0 then equation (11) becomes: ( ) ( ) 2 2 2 0 2 1 λ ζ λ + - = B B . ( 12
) or ( ) ( ) 2 2 2 0 2 1 λ ζ λ + - = iKB F . ( 13
)
If we define 0 B B = θ then equation (12) becomes:

( ) ( ) 2 2 2 2 1 1 λ ζ λ θ + - = . ( 14
)
Differentiating equation ( 14) with respect to λ and setting the derivative equal to zero, we obtain:

2 2 1 ζ λ - = . ( 15
) or 2 2 2n p - = ω . Q . 2 2 2 p n << ∴ m K i = ω . (16)
Equation ( 16) is the requirement of occurring resonance for the device.

Based on damping theory the damping coefficient is [START_REF] Zhang | Mechanical Vibration[END_REF]:

B F e ω π γ ′ = 4 . ( 17
iKm B F ω π ζ ′ = 2 . ( 18
)
Substituting equation ( 15)and(18) into equation (13) yields:

( ) 2 2 2 0 4 ) ( 4 F Km i B mB F F ′ - ′ = ω π ω π . ( 19
)
Equation ( 19) is used to evaluate the amplitude of the exciting force. Equation

(2), ( 16) and ( 19) is applied to evaluate the magnitude of the exciting force.

Analysis of Damping

The forces exerted on the rectangular screen at the plane perpendicular to the vibrating direction are(Fig. 3): composite pressure force of the two 5-dampers, ;

composite pressure force of the two 6-dampers, ; reacting forces of the frame, and ; and component of gravitational force of the rectangular screen that includes the tuber and soil on it, The composite pressure force of the two 5-dampers is determined as:

1 N 2 N 1 N ′ 2 N ′ β sin mg . 2 N β sin mg 1 N ' 2 N ' 1 N rectangular screen
1 1 1 2 x k N = . ( 20
)
Where is stiffness of the spring in the 5-damper ( ), is the deformed

1 k 1 - m N 1 x
magnitude of the spring in the 5-damper ( ).

m

The composite pressure force of the two 6-dampers is determined as:

2 2 2 2 x k N = . ( 21
)
Where is stiffness of the spring in the 6-damper ( ), is the deformed magnitude of the spring in the 6-damper ( ).

2 k 1 - m N 2 x m
The reacting forces may be expressed as:

1 1 1 2 x k N = ′ . ( 22
)
β sin 2 2 mg N N + = ′ . ( 23
)
The damping may also be evaluated by:

2 2 1 1 N N N N F ′ + + ′ + = ′ μ μ μ μ . ( 24
)
Where μ is friction coefficient between the 5-damper or the 6-damper and the rectangular screen . It is also the friction coefficient between the slideways of the frame and the rectangular screen .

Substituting equation ( 20

g m F x k x k - ′ = + . ( 26
)
Equation ( 26) is applied to evaluate the deformed magnitude of the damping springs in the four dampers.

Results and Discussion

A motion simulation experiment for the rectangular screen on a potato harvester has been done with the Matlab/Simulink module. The parameters of the

Fig. 1 .

 1 Fig. 1. A model of a soil-tuber separation device: (a) front view; (b) K direction view except the linear motor .1-linear motor; 2-rectangular screen; 3-supporting springs; 4-frame; 5-dampers; 6 -dampers. arrow denotes direction of travel .

.

Fig. 2 .

 2 Fig. 2. Forces acting on a rectangular screen in vibrating direction.

Fig. 3 .

 3 Fig. 3. Forces acting on a rectangular screen at the plane perpendicular to the vibrating direction.

  ),(21),(22),and(23) into equation (24) and rearranging terms gives:

 2), ( 16) and ( 19). The simulation result was as shown in Fig. 4. We could see that the motion trail of the rectangular screen is correct. In practice, the stiffness and the deformed magnitude of the springs in the dampers could be adjusted. Therefore,the appropriate damping (2) Forces acting on the rectangular screen are found. A differential equation for describing the forces is set up. An equation for evaluating exciting force of the linear motor has been derived from the solution to the differential equation. An equation of the factors affecting on the device resonance has also been obtained.

Additionally, an equation for evaluating the deformed magnitude of the damping springs in the dampers has been derived.

(3) A motion simulation experiment for the rectangular screen on a potato harvester has been done with the Matlab/Simulink module. The experiment expressed that the method is feasible.