{'fraction_non_alphanumeric': 0.032404028608962195, 'fraction_numerical': 0.009779594219821924, 'mean_word_length': 4.554249425753277, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Introducing Randomized High Order Fuzzy Cognitive Maps as Reservoir Computing Models: A Case Study in Solar Energy and Load Forecasting

Fuzzy Cognitive Maps (FCMs) have emerged as an interpretable signed weighted digraph method consisting of nodes (concepts) and weights which represent the dependencies among the concepts. Although FCMs have attained considerable achievements in various time series prediction applications, designing an FCM model with time-efficient training method is still an open challenge. Thus, this paper introduces a novel univariate time series forecasting technique, which is composed of a group of randomized high order FCM models labeled R-HFCM. The novelty of the proposed R-HFCM model is relevant to merging the concepts of FCM and Echo State Network (ESN) as an efficient and particular family of Reservoir Computing (RC) models, where the least squares algorithm is applied to train the model. From another perspective, the structure of R-HFCM consists of the input layer, reservoir layer, and output layer in which only the output layer is trainable while the weights of each sub-reservoir components are selected randomly and keep constant during the training process. As case studies, this model considers solar energy forecasting with public data for Brazilian solar stations as well as Malaysia dataset, which includes hourly electric load and temperature data of the power supply company of the city of Johor in Malaysia. The experiment also includes the effect of the map size, activation function, the presence of bias and the size of the reservoir on the accuracy of R-HFCM method. The obtained results confirm the outperformance of the proposed R-HFCM model in comparison to the other methods. This study provides evidence that FCM can be a new way to implement a reservoir of dynamics in time series modelling.


Introduction
In the past decades, time series analyzing and forecasting, as a hot research topic, has attracted the attention of many researchers to make better decisions across various disciplines specially with the advent of some recent technologies such as loT, IoE or Big Data. Despite the proposal of plenty of various forecasting procedures in the literature, the accurate prediction of the future values is still an open challenge due to the existence of non-linearity and uncertainty in real time events. Statistical time series forecasting methods are limited due to some characteristics such as time consuming, lack of scalability and explainability, inability to deal with uncertainty and complex problems in the real world [1]. Accordingly, fuzzy time series (FTS) was presented by Song and Chissom in [2,3,4] to handle the mentioned problems. Numerous FTS forecasting models have been introduced in the literature in many applications due to some important features including simplicity, explainability, flexibility, updatability, readability, versatility, and high accuracy of forecasting [1,5].
Fuzzy Cognitive Map (FCM) was proposed by Kosko in [6] as weighted knowledge-based model that is used to do the task of forecasting by extracting knowledge in FTS models. FCM is a kind of interpretable recurrent neural network that combines the concepts of fuzzy logic and neural networks which is composed of nodes (representing concepts) and signed directed relations (causal relations) between a couple of concepts. Noteworthy that FCMs as qualitative soft computing techniques can be used to represent the dynamic behaviour of complex systems with high ability to dealing with uncertainties [7]. Therefore a wide range of FCM-based time series forecasting methods have been developed in the literature [8,9,10,11,12,13,14,15,16,17].
Based on the literature, GA and PSO have been used widely to learn weight matrices in FCM-based time series forecasting models using evolutionary-based learning methods. For example, [25] proposed a univariate forecasting method using FCM and C-means clustering plus moving window technique, while PSO algorithm was used to learn weight matrices. In [11], the authors proposed univariate time series forecasting using FCM and C-means clustering in which real coded genetic algorithm (RCGA) was employed to extract weight matrix. [20] introduces a time series forecasting model based on the synergy of High-Order FCM (HFCM) and fuzzy C-means clustering in which PSO is applied to optimize the weights of HFCM. In [8] a two-level forecasting technique was proposed to carry out forecasting in both numerical and linguistic terms using RCGA learning method to train weight matrices. A double-phased approach was introduced in [24] combining granular FCM and fuzzy C-means clustering which is trained through GA algorithm. The authors in [38] proposed univariate HFCM-FTS model using GA to train weight matrices. Intuitionistic fuzzy grey cognitive maps for forecasting interval-valued time series learned via Differential Evolution (DE) was proposed in [39].
Since the population-based learning methods are much time consuming and computationally expensive due to the large number of learning parameters and large number of learning processes, finding and designing fast and robust learning strategies is considered the main challenges in this area. Therefore, evolutionary learning has been replaced with other techniques in some references. For instance, ridge regression in [13,17,40], Bayesian ridge regression in [21], Moore-Penrose inverse in [16]. A rapid and robust learning method with maximum entropy was proposed in [41] to learn large scale FCMs composed of least-squares and maximum entropy terms. The robustness of the well-learned FCM is guaranteed by Least-squares term and the maximum entropy term regularizes the distribution of the weights of the well-learned FCM. Recently, a new time series modeling based on least square FCM termed as LSFCM was introduced in [42]. In this method Fuzzy c-means clustering is exploited to construct FCM concepts and the weight matrices are adopted from the given historical observation of time series using the least square method which is much faster than other population-based methods.
The crux objective of our research study is to develop a novel univariate FCM-based time series forecasting model focusing mainly on FCM training to promote the accuracy and efficacy of time series methods applying FCM. Thus, the concepts of High Order FCM (HFCM) and Echo State Network (ESN) reservoir computing are merged to construct a novel forecasting technique using the least square training algorithm termed as Randomized HFCM (R-HFCM). More vividly, this model consists of three layers: input layer, reservoir and output layer. The reservoir layer consists of a group of HFCM-FTS models proposed in [38] but in this case the weights are randomly initialized such that the Echo State Property (ESP) condition in ESN reservoir computing is satisfied. Then, each sub-reservoir generates its output independently after feeding input into each sub-reservoir separately. Finally, the least squares minimization technique is applied to train the output layer and generate the final predicted value. As the computational experiments reveal, in addition to improving the accuracy, the proposed model is much faster than HFCM-FTS method. Furthermore, the results are promising in comparison to other state-of-the-art techniques in the literature.
The remainder of this paper is organized as follows: Section 2 presents a brief description of the Fuzzy Cognitive Map (FCM) and Echo State Network (ESN); Section 3 introduces the proposed method in details; Section 4 presents the experimental results and discussion and finally the paper conclusion and some possibilities of future works are drawn in Section 5.

Fuzzy Cognitive Maps
Fuzzy Cognitive Map (FCM) was proposed by [6] as qualitative and causal model which combines fuzzy logic and CM to present uncertainties and complex characteristics of the systems. Graphically speaking, FCM is powerful interpretable knowledge-based model which is composed of concepts (nodes) and the causal connection among the pair of concepts . The directed signed arrows among the concepts is known as weights which reflect the effect of one node on another one. Therefore, FCM is based on cause and effect causality. If the node influences on others is cause and while influenced by others is an effect.
Each FCM is identified by four elements which is defined by 4-tuple (C, W, a, f ), where C = [c 1 , . . . , c n ] is the set of n concepts, which are the variables (and the nodes of the graph) that compose the system. The state values (activation degree) of these concepts at any time t is represented as follows: a = (a 1 , . . . , a n ) (1) where a i ∈ [0, 1], i = 1, 2, ..., n. The relations between the nodes are described as an n × n weight matrices as follows: where w ij ∈ [−1, 1] (i, j = 1, 2, . . . , n) indicates the relationship between the source and the target nodes which can be positive, negative or zero. If w ij > 0, a decrease/increase in value of node c i makes a decrease/increase in value of node c j . If w ij < 0, a decrease/increase in value of node c i makes an increase/decrease in value of node c j . If w ij = 0, there exists no relation between nodes c i and c j . Figure 1 exhibits the simple example of FCM structure including 5 nodes. The final element is activation function (f ) used in Kosko's updating rule to preserve the activation degree of each concept inside the predefined interval. Bivalent, trivalent, hyperbolic tangent and the sigmoid are the most common activation functions in the literature [23]. According to the Kosko's activation rule, the activation degree of each node at t + 1-th iteration is updated considering the weight matrix and activation degree of all connected nodes at the t-th iteration. In other words, the weights are used to transfer the current activation state of each concept to the next. Hence, the following formula describes the dynamics of FCM.
(3) where w ij highlights the value of the causal relationship between concepts c i and c j , whereas a t i represents the state value of concept c i at time step t.
Although the aforementioned updating rule has been utilized widely in many FCM-based applications, other modifications have been proposed in the literature with/without self-connection and with/without considering memory [43,23,44,45,46].
Dynamically speaking, the equation (3) expresses only the first order dynamics of the FCM. It means that the value of each concept at time t + 1 only relies on the activation level of all concepts at time t. But an accurate modeling of FCM will not be obtained only by considering the current activation value of the concepts and ignoring the past values. Accordingly, High-Order FCMs (HFCMs) were proposed to improve the performance of FCMs and cover the mentioned limitation. Thus, the equation (3) is modified to describe the dynamic behavior of complex systems more accurately. The below equation describes the k-order FCM: [13,20] where w 0 i stands for the bias term and w Ω ij is the casual relation originating from c i and pointing to c j at time step t − Ω + 1. Based on this equation, the activation level of i-th node at the moment t + 1 depends on the activation degree of all concepts at {t, t − 1, .., t − Ω + 1} moments, not only the activation states of the concepts at time t, during the iterative process.
Equation (4) can be rewritten in matrix form as:

Echo State Network (ESN)
Echo state network (ESN) was proposed by Jaeger [47] as one of the improved Recurrent Neural Network (RNN). Standard RNN has been considered as a powerful tool to simulate complex dynamic systems, but training algorithms of RNNs involve some downsides including: relatively high computational training costs and potentially slow convergence, local minima of the error function and vanishing of the gradients. With the goal of handling these issues, a new approach in RNN training called Reservoir Computing (RC) was proposed to study of initialization conditions and stability instead of focusing on the training algorithm. ESN as one of the popular RC techniques, has a very simple and fast training process and it has abundant nonlinear echo states and short-term memory because of high-dimensional projection and highly sparse connectivity of neurons in the reservoir, which are useful for modeling dynamical systems [48]. The structure of ESN consists of three basic components including input layer, a large recurrent hidden layer (called the dynamical reservoir as untrained large, sparsely connected nonlinear layer) with fixed sparse hidden-to-hidden connections and an output layer. Despite traditional RNNs, the dynamical reservoir or hidden layer in ESNs is untrainable and only the output weights are trained which is the most excellent feature of ESNs. Figure 2 highlights the structure of traditional ESN with M input units, N reservoir units and L output (readout) neurons, while u(t) = (u 1 (t), . . . , u M (t)) T , x(t) = (x 1 (t), . . . , x N (t)) T , and y(t) = (y 1 (t), . . . , y L (t)) T determine their activation at each time step. The output of dynamic reservoir neurons and the output of ESN at time step t + 1 are updated according to the following equations respectively.
where the N × M input weight matrix W in represents the relation among input units to the reservoir neurons; the internal connection among neurons inside the reservoir layer is represented via N × N weight matrix W res ; the output weight matrix is represented via W out with L × (M + N ) dimension. Finally, an N × L weight matrix represents the connection projected back from readout neurons to reservoir units as W fb or W back . Interestingly, both internal and input weights are initialized randomly with no changes during the training and testing process. On the other side, a simple linear regression is applied to train readout layer. Noteworthy that an ESN is trained through a supervised learning method in two stages. At the first level, the M dimensional inputs are mapped into a high dimensional reservoir state to reach the echo states x. Then a simple regression method is exploited to train the output weights.
In addition to the mentioned unique features of ESN, regarding the fixed non trainable reservoir and just learning readout layer, the ESN can also act as a kernel in kernel-based learning technique to catch the input dynamics by application of a random high-dimensional projection method which individuates the ESN from other RNNs [48]. In a nutshell, the learning of ESN is comparable with other RNNs because it is simple and fast. Besides, it does not trap into local minimum which enables ESN for modeling dynamics of time series. Due to the simple method and high learning efficiency, ESN has been successfully applied to many fields, such as time series prediction tasks [49,50,51], dynamic pattern classification [52,53], speech recognition [54,55], and so on. In other words, a very easy training process, high-dimensional projection as well as highly sparse connectivity of neurons in the reservoir enable ESN to modeling dynamic systems.

Initialization and Hyper-parameters
An ESN is constructed according to the significant hyper-parameters including the reservoir size N , the input scale IS, the spectral radius , and sparsity γ. The input scale is employed to initialize the input to hidden weight matrix W in such that each member of W in obeys the uniform distribution in [−IS, IS]. γ determines the proportion of non-zero elements in (W res ) and is the spectral radius of W res , which must be set smaller than 1 [48].
As mentioned, in the training process of ESN, just the readout weight matrices are trainable while the W res is selected randomly with a uniform distribution symmetric around the zero value before training execution to provide the requirements of echo state condition [47,48]. The initialization of W res is done through the following procedure to assure that the maximum absolute eigenvalue or spectral radius is less than one. In other words, firstly an internal randomly weight matrix W rand will be generated. Then rescaled to meet the Echo State Property (ESP) condition as below: where ρ max (W rand ) denotes the maximum eigenvalue of matrix W rand when its elements are generated randomly. To guarantee the stability of ESN, must be set smaller than 1. This is a necessary condition of ESN stability [47]. The final step is to determine W out which is unknown and alterable. W out is calculated through the following learning equation: where (X T X) −1 denotes the inversion of square matrix X T X. Noteworthy, ESNs with leaky integrator neurons [56], φ-ESN in [57], ESNs with circle reservoir topology [58] and different types of deep ESNs in [59,48,60,61] represent powerful variations and advances in ESN in the literature.

Proposed R-HFCM method
As was mentioned earlier, the large number of processes, time-consuming and learning parameters adjustment are counted as the major deficiencies of the population-based FCM learning methods. Thus, to rectify the above-mentioned limitations, designing a faster training model seems vital. Due to that, some publications focused on other new techniques which are much faster. For instance, the authors in [42] have developed a new time series modeling based on least squares FCM termed as LSFCM that is much faster rather than multi-iteration stochastic searching. Accordingly, there is no exception for the proposed HFCM-FTS model in [38] which has been trained by GA. Therefore, this part introduces novel methodology to overcome the mentioned fundamental problem. This model is composed of a group of randomized HFCM-FTS models termed as Random HFCM (or R-HFCM) model to predict univariate time series. The innovation of the proposed R-HFCM model is integrating the concepts of FCM and ESN reservoir computing method exploiting the least-squares regression algorithm as the learning strategy which is detailed in the following.
From a structural perspective, this model consists of three layers as illustrated in the Figure 3 including the input layer, the intermediate or reservoir layer and the output layer. The reservoir is composed of a group of random HFCM-FTS models that are fed only by the external input time series to provide inputs for the output layer. Then the obtained outputs from each sub-reservoir are considered as the units of the output layer that are used to generate the final predicted value. As deduced from Figure 3, the reservoir structure consists of some numbers of sub-reservoirs (N SR ) so that there are no relations among them. It means that the proposed model does not have a hierarchical stacking structure similar to the deep models. Since the structure of each sub-reservoir in the proposed R-HFCM model exactly is designed the same as the HFCM-FTS model, the R-HFCM model is mainly concentrated on training the output layer using least squares regression. From another perspective, it can be said that the proposed R-HFCM is a kind of ESN in which only the output layer is trainable while the reservoir parameters are initialized randomly and remain unalterable during the training process. On the other hand, the weights of each sub-reservoir in the reservoir layer are chosen randomly to meet the ESP condition in ESN and then each sub-reservoir generates its output independently. It means that in our model the reservoir weights are initialized based on the ESP condition in ESN reservoir computing, then the least squares minimization algorithm is applied to the output units to generate the final predicted value.
The proposed R-HFCM method is divided into two procedures: the Training procedure and the Forecasting procedure. Therefore, the Training and the Forecasting procedures are detailed in the rest of this section.

Training Procedure
The main objective of the training procedure is to generate the linguistic variable C with k concepts, to find the best least square coefficients to fit the data, given a crisp training set Y and the activation function f informed by the user. The steps of the method are listed below:

Partitioning:
Since R-HFCM is composed of a group of HFCM-FTS models, the uniform scheme is applied for partitioning process in both R-HFCM and HFCM-FTS models but in this case for each sub-reservoir. More precisely, the Universe of Discourse (U), for each sub-reservoir, is partitioned into k even length and overlapped intervals. Then, a fuzzy set C i (i.e. concept) is defined with a membership function µ Ci . In each sub-reservoir, the group of the k concepts form the linguistic variable C, such that C i ∈ C, ∀i = 1, . . . , k. Here the grid partitioning is used to generate versions of the FCM with k = {5, 10, 20}, where the number of concepts is equal to the number of partitions and triangular membership function µ Ci . Figure 4 shows a simple example of the FCM structure used in this model. In Figure 4-A the universe of discourse is partitioned to generate seven concepts, which is also the number of fuzzy sets because in our model the number of concepts and fuzzy sets is considered the same. As shown in Figure 4-C, FCM is a collection of nodes (concepts) and causal interactions among these concepts which is represented via the weights matrix. According to Figure 4-B, this weight matrix is a square connection matrix which is defined randomly through the following step.

Least Squares coefficients determination:
Each matrix W t , for t = 1, . . . , Ω, is a k × k matrix where w t ij ∈ R is the weight between the concepts C i and C j at the time lag t, and Ω is the order of the model. As noted earlier, the innovation in the proposed model is that the weight matrices of each HFCM-FTS are randomly chosen from a uniform distribution over interval [-1,1]. Then, they are scaled according to the ESN reservoir computing to preserve the ESP condition as described by the following formula: where ρ max (W rand ) is the maximum eigenvalue of W rand and ∈ (0, 1) is the desired spectral value (scaling parameter). Also, since the members of the weight matrix in the FCM must be in the range [-1,1], the value of in our model is equal to 0.5 to satisfy the corresponding condition. Hence, each reservoir generates its own output through the Forecasting procedure, given the training sample as well as the initialized weights. Finally, the least squares regression method is applied to train the last layer. Through this way, a linear model is solved and the optimum least-squares coefficients are obtained by minimizing the error function. Figure 5 shows a simple structure of R-HFCM model considering N SR = 2 while in each sub-reservoir the order is Ω = 2. Thus there are two weight matrices for times t and t − 1 for each randomized HFCM-FTS as indicated in the figure. It is worth noting that the model contains the bias weights while they are discarded from the Figure 5 just because of ease of notation.

Forecasting Procedure
The main goal of this stage is to compute the predicted crisp valuesŷ(t + 1) of each sub-reservoir as well as the final predicted valueŷ f (t + 1) by given the linguistic variable C, the weight matrices W t , the activation function f and a crisp input Y (t). The steps of this stage are listed below:

Fuzzification:
In each sub-reservoir, given the crisp input sample Y with size T , each instance y(t) ∈ Y, t = 1..T , is transformed into an activation vector a(t) such that a i (t) = µ Ci (y(t)), ∀C i ∈ C, that is, each value a i (t) ∈ a(t) corresponds to the membership degree of y(t) to the concept C i .

Activation:
For each HFCM-FTS, the state value of each concept in time t + 1 can be defined by the following formula: Since this model considers the presence of the bias term, w 0 represents bias in the above equation.
It should be noted that the initialization of the bias parameters and weight matrices are the same, such that the ESP condition in ESN reservoir computing is satisfied. It means that they are randomly chosen from a uniform distribution over the interval [-1,1], then rescaled according to the equation (10) but in this case, the weight matrices are replaced by the bias ones.

Defuzzification:
After calculating the activation level of each concept for each sub-reservoir, in this step, the defuzzification is carried out to produce the output related to each of them. Therefore, the forecast values for each sub-reservoir at time t + 1 in numeric terms can be calculated via the below equation: where a i (t + 1) is the activation calculated from the previous step for each concept at time t + 1 and mp i is the center of each concept C i . Finally, the predicted value at time t + 1 is computed through the linear combination of the obtained outputs of all sub-reservoirs and the least squares coefficients. The final predicted value is described as following: For instance, as highlighted in Figure 5, the final output is estimated via the below equation when N SR = 2 :ŷ f (t + 1) = λ 0 + λ 1 ·ŷ 1 (t + 1) + λ 2 ·ŷ 2 (t + 1) where {λ 0 , λ 1 , λ 2 } denote the obtained least square coefficients from training procedure andŷ 1 (t + 1) andŷ 2 (t + 1) are the generated outputs for the first and second sub-reservoirs respectively according to the equation (12).
As Figure 5 recounts, the number of least squares coefficients directly depends on the number of subreservoirs N SR . To be more precise, the number of least squares coefficients in this model is equal to N SR + 1. Thereby, for N SR = n, the final predicted value is described by the following formula using n + 1 least-squares coefficients. y f (t + 1) = λ 0 + λ 1 ·ŷ 1 (t + 1) + λ 2 ·ŷ 2 (t + 1) + ... + λ n ·ŷ n (t + 1)

Dataset
In order to evaluate the effectiveness of the proposed R-HFCM model, in this section, two different datasets including the SONDA dataset and the Malaysian dataset are employed.

Variable
Type Description DateTime Time Stamp yyyy-MM-dd HH:MM glo avg Real Global average solar radiation ws 10m Real Wind speed in meters by second (m/s)  Figure 6 shows the weekly and yearly samples of the referred solar radiance time series. Also, Figure  7 highlights autocorrelation and partial autocorrelation plots for solar radiance time series. In this experiment, as shown in Figure 8 The performance of the proposed R-HFCM model is influenced by some factors. The number of concepts (partitions), the type of activation function, bias terms and the size of reservoir (the number of sub-reservoirs N SR ) are the significant parameters affecting the accuracy of the proposed method. Thus, in the following sections, we assess the performance of the proposed model for SONDA and Malaysia datasets respectively with respect to the mentioned influential parameters and compare the results with some recent state-of-theart time series forecasting methods.

SONDA Case Study (A) Results and Discussion
In the first scenario, 8000 samples of the SONDA dataset have been used as can be seen from Figure 8. A sliding window of 2, 000 samples in the cross validation method has been used. 80% of the window  Table 3, different activation functions were used to design R-HFCM model. Table 4     On the other side, with a constant value of k, the error rate follows a downward trend as the number of N SR increases. Thus, the model is highly sensitive to the selection of both k and N SR parameters and the optimal model is designed by increasing N SR and inserting 5 concepts in each sub-reservoir. Also, as depicted in the figures, the effect of activation function over the model accuracy is less than the other factors. This means that there is no significant difference in the accuracy of the model when using different activation functions, especially in the case of using higher N SR and lower k. Therefore, as presented in Figure 15, the model performance is less sensitive against the choice of the activation function. In this figure, each box contains six values in terms of average RMSE and U for different values of N SR with regards to the specified k number and a certain type of activation function.

(B) Comparison with other Fuzzy Time Series methods
In this part, the proposed method is tested against other recent fuzzy time series forecasting methods available in the literature. In particular, the proposed model is tested against HFCM-FTS , HOFTS, WHOFTS and PWFTS methods. Bear in your mind that since the weights and bias are randomly chosen according to the ESN reservoir computing to meet the ESP condition, the R-HFCM results are obtained after 20 independent runs. Table 5 shows the data extracted from Table 4 on the basis of the best (minimum) average RMSE and U. Accordingly, Table 5 compares the best results obtained for each activation function with respect to the different number of concepts. More detailed look at the results reveal that with the same map size k, the best performance can achieve with the same or different size of reservoirs. For example, when k = 5 , the reservoir size of the models with ReLU and Softplus activation functions are respectively greater and equal to the size of the reservoir of the model using sigmoid. Although there exist slight differences among the model accuracy using different activation functions, irrespective of the number of concepts, the minimal error is achieved by considering Softplus as the model activation function. More precisely, selecting this activation function leads to the best model performance with unequal size of reservoirs such that for k = {5, 10, 20} concepts, the optimum models are designed utilizing N SR = {20, 30, 40} sub-reservoirs, respectively. Therefore, the obtained results using Softplus are exploited to compare the suggested model with other recently presented methods in the literature. Accordingly, Table 6 compares the results of the proposed methodology with the aforementioned methods with respect to the effective parameters. The results indicate the significant supremacy of the proposed R-HFCM over HFCM-FTS model. Actually, in comparison to the HFCM-FTS model, the proposed model is much faster and more accurate with fewer parameters. In other words, since the weights and bias matrices are determined randomly in the proposed model, the number of parameters only depends on the number of N SR . Thus, the least squares regression algorithm in the proposed model will solve a linear problem with fewer variables compared with using GA in HFCM-FTS model. For instance, with 10 concepts, the HFCM-FTS parameters are more than 7 times that of the proposed R-HFCM model which can be considered as an effective factor on the training time of the model. Interestingly that by increasing the number of concepts from 10 to 20, only 10 parameters have been added to the proposed model, while this number will be almost quadrupled in the HFCM-FTS model. Additionally, adjustment of the learning parameters in GA is considered as another limitation compared with least squares method. Thereby, the training time of the model will be significantly affected by the mentioned factors.

Malaysia Case Study
In the second scenario, in order to appraise the validity of the proposed model, we consider the Malaysia data set, which includes hourly electric load and temperature data. 8760 samples of load and temperature are selected as the model inputs as shown in Figures 11 and 12 respectively. Since the model has been designed to predict univariate time series, hourly electric load and temperature data are fed into the model separately. Thus, the rest of this section examines the obtained results using hourly electric load and temperature datasets, respectively.

(A) Hourly electric load dataset
Similar to the first case study, a sliding window of 2, 000 samples in the cross validation method has been used but in this scenario 90% of the window for training and 10% for test. The experiment has been executed with k = {5, 10, 20} concepts, N SR = {2, 5, 10, 20, 30, 40} sub-reservoirs, Softplus activation function and by considering bias terms to assess the model performance. As illustrated in Table 7, RMSE, MAPE and U are considered as forecasting accuracy metrics to assess the model performance. So it shows the average of each metric for the 20 independent runs given different numbers of concepts as well as reservoirs of different sizes.
As the table summarizes, like SONDA case, the accuracy of the model has direct and indirect relation with N SR and k respectively. The least value of the average error is generated with higher values of N SR for model with 5, 10 and 20 concepts, respectively. Also, Figure 16 highlights clearly the relation between the model accuracy with N SR and k in terms of the average MAPE and U by using Softplus activation function. More specifically, the error rate is decreasing from left to right. In simple terms, the higher value of N SR , the more precise the forecasting model becomes. For instance, in all cases, the performance of the model is more accurate using N SR = 40 compared to the model with N SR = 5. Thus, with a fixed number of concepts, the higher the N RS , the lower the forecasting error. In reverse, the lower the number of partitions, the better the model accuracy. It means that with a constant N RS , increasing the number of k leads to generate the model with less accuracy. Therefore, the lowest error rate is obtained with the minimum value of k and with increasing the size of the reservoir. Noteworthy, for k = {5, 10} and N SR ≥ 20, the performance of the models are different slightly. For instance, for k = 5 the model performances with N SR = {30, 40} are very close. (B) Hourly temperature dataset In this section, the codes are implemented to evaluate the model efficacy, while hourly temperature data is fed into the model as input instead of hourly load. In this case, like hourly load data, 2, 000 samples in the cross validation method have been used, with the same segmentation for training and testing. Also, this experiment is carried out with the same conditions, i.e., k = {5, 10, 20}, N RS = {2, 5, 10, 20, 30, 40}, Softplus activation function and the existence of bias terms. According to Table 8, RMSE, MAPE and U are considered as the accuracy metrics to evaluate the performance of the proposed method after 20 independent runs. The table records the average values of each accuracy metric applying the defined values of k and N SR . As can be seen from Table 8, similar to the load case study, the relation between the model accuracy with k and N SR follows the uniform pattern. That is, the fewer k number and the more N SR number, the more accurate the model. Therefore, with a constant concept number k, the accuracy is improved by increasing the number of N SR . For instance, assume that k = 5, the model with N SR = 30 is much more accurate than the model with N SR = 2. On the other side, with a certain number of N SR , boosting the number of concepts decays the model accuracy performance. For example assuming N SR = 20, the accuracy of the model is lessened by raising the value of k from 5 to 20.
These relations between the model accuracy with the effective elements on the model accuracy can also be deduced from Figure 17. It highlights that for any number of concepts the average errors are reduced by increasing the size of reservoir. Although the higher accuracy goes to the model with minimum k,

RMSE
Ave.  In addition, Figure 18 shows the box-plot that reflects the sensitivity of the method accuracy to the effective parameters including the map size and reservoir size for both load and temperature dataset. Each box consists of six values including the average error in terms of MAPE or U considering N SR = {2, 5, 10, 20, 30, 40} with respect to the defined number of concepts. As the figure recounts, increasing the number of k deteriorates the accuracy. Overall, by considering the combination of the number of concepts and sub-reservoirs, it is clear that the best performance is achieved for the least values of k and the higher number of N SR . Accordingly, Table 9 summarizes the best model performance for both load and temperature data in terms of the best average errors. As the table confirms, by increasing the number of concepts, the model accuracy is declined. Also, for different number of concepts, the best results are obtained with the same or different size of reservoirs.

(C) Comparison With Other Methods
In this section, the proposed R-HFCM model is tested against other investigated models in [62] under the assumption that the goal is to predict the hourly electric load. Accordingly, Table 10 shows the comparison among the obtained results from the proposed R-HFCM method and other models including FTS-CNN, LSTM, PWFTS, and SARIMA in terms of average MAPE. R-HFCM is run with k = 5, 10, 20 concepts. As the table suggested, the performance of the R-HFCM (k = 5) and LSTM are very close in terms of accuracy. It is interesting that, although the results indicate the superior performance of the FTS-CNN model over other methods, the results of the proposed R-HFCM model are still promising. This is because unlike FTS-CNN and LSTM methods, R-HFCM is a univariate model that utilizes only one variable (load) to predict hourly electric load. Also, R-HFCM model is a shallow model with lower parameters and cheap computational cost in comparison to the FTS-CNN and LSTM methods, which are deep learning models. Briefly, the proposed R-HFCM model which is equipped with a fast and nonrepetitive learning algorithm is competitive with other FTS methods and some other state-of-the-art methods.
In addition, in order to compare the performance of the proposed R-HFCM model, Figure 19 represents the Friedman Aligned Ranks of the methods. The test statistic for these results is Q = 6.0 where the p-value is 0.42319. For this statistic value, the H 0 is accepted at α = 0.05 confidence level which indicates that there is no difference between the means of the competitor models. In other words, all the methods are statistically equivalent. Therefore, the R-HFCM (k = 5) method performed satisfactorily when compared with the standard methods in the literature.

Conclusion and Future Works
This paper provides a new randomized-based method to predict univariate time series using FCM. It is conducted based on the introduction of a novel FCM forecasting model in the literature termed as R-HFCM. R-HFCM model is a group of random HFCM-FTS models integrating the concepts of FCM and ESN. To be more clear, R-HFCM is composed of three layers employing random HFCM-FTS models as the components of the reservoir unit called sub-reservoirs. Since the structure of each sub-reservoir is exactly the same as HFCM-FTS, the main focus of this model is mainly on training the output layer using least squares regression. Thus, the weights of each sub-reservoir in the internal layer are initialized randomly so that the ESP condition in ESN is met and produces the output from each sub-reservoir. Then least squares regression method is applied to the output units to generate the final predicted value. Also, this method considers the effect of some elements on the model accuracy including the number of concepts, the number of sub-reservoirs, activation function and the presence of bias term. Based on the results, the more sub-reservoirs and the fewer concepts, the more accurate the model. That is the best result is achieved by increasing N SR and considering the minimum number of concepts in each sub-reservoir. Meanwhile, the activation function has the least effect on the accuracy of the proposed method compared with N SR and k.
It is considerable that R-HFCM offers a mixture of the advantages of high speed and accuracy over the HFCM-FTS model while the relation between the accuracy and the number of concepts is similar in both models. This means that as the number of concepts increases, the accuracy of the prediction, the interpretability and readability of the model deteriorates. Moreover, the number of parameters in R-HFCM directly depends on N SR . Thereby, the least square method will solve a linear problem with fewer variables than GA, which can be counted as an effective element on the training speed of the model.
In the proposed R-HFCM model, the reservoir components are fed only through the external input to generate the input for the output layer. Then, the obtained results from each sub-reservoir are exploited to make the final predicted value. In this way, the proposed R-HFCM model has a shallow structure with no relations among the sub-reservoirs. The generated output from each layer of reservoir is only fed to the output layer of the model. Accordingly, the future challenge may involve extending the proposed model from shallow to deep one. Furthermore, as the model has been designed to predict univariate time series, predicting high-dimensional multivariate time series can be considered as the other possibility for future work.