extfieldsofstudy:
Computer Science
Medicine

text:
Method for Handling Massive IoT Traffic in 5G Networks

The ever-growing Internet of Things (IoT) data traffic is one of the primary research focuses of future mobile networks. 3rd Generation Partnership Project (3GPP) standards like Long Term Evolution-Advanced (LTE-A) have been designed for broadband services. However, IoT devices are mainly based on narrowband applications. Standards like LTE-A might not provide efficient spectrum utilization when serving IoT applications. The aggregation of IoT data at an intermediate node before transmission can answer the issues of spectral efficiency. The objective of this work is to utilize the low cost 3GPP fixed, inband, layer-3 Relay Node (RN) for integrating IoT traffic into 5G network by multiplexing data packets at the RN before transmission to the Base Station (BS) in the form of large multiplexed packets. Frequency resource blocks can be shared among several devices with this method. An analytical model for this scheme, developed as an r-stage Coxian process, determines the radio resource utilization and system gain achieved. The model is validated by comparing the obtained results with simulation results.

Introduction
Internet of Things (IoT) communication has gained significant popularity in recent years. Existing IoT communication is reliant on present-day mobile networks. The current networks are meeting the requirements of existing IoT devices sufficiently to date. However, the massive growth of IoT data traffic means that the existing mobile communication systems (3G, 4G, etc.) might not be capable of coping with substantial future data traffic. IoT traffic is usually dissimilar to regular traffic such as conversational video or file transfer. IoT data are mostly generated by a large number of IoT devices in the form of small packets. IoT communication is mostly based on narrowband applications. Conversely, the human-generated normal data traffic appears from small number of mobile phones in the shape of large sized packets. Meanwhile, 5G is now being considered as a strong candidate for serving the massive future IoT traffic. Rollouts of 5G networks are expected in the near future.
The 3rd Generation Partnership Project (3GPP) standard, Long Term Evolution-Advanced (LTE-A) has been developed for broadband applications. With narrowband applications, LTE-A is not capable of achieving efficiency in terms of bandwidth usage and cost. Incorporation of the narrowband IoT data traffic into LTE-A system requires revision in the design of such networks. Otherwise, considerable degradation of the overall network performance can be faced. Therefore, achieving for modelling of multiplexing buffer. Authors of [23] present a mathematical model for evaluation of AAL2 type multiplexer's performance parameters like expected sojourn time, mean queue length, and the probability of delay violation. Authors of [24] analyze the AAL2 multiplexer by modelling the multiplexer departure process as an r-stage Coxian process. The authors of [25] present the performance evaluation of voice traffic in mobile radio networks with Coxian distributed channel holding time. The authors of [26] propose a method of performance analysis for a continuous time Markovian arrival process at the ATM Adaptation Layer 2 (AAL2), where the Common Part Sublayer (CPS) packets get multiplexed at the AAL2 buffer before being sent to transmission buffer. Both the buffers are modelled as queues.

Air Interface of 5G
The primary motivation behind the development of new standards for 5G mobile systems is innovations in support of IoT traffic. 5G applications for IoT are expected to increase market space and open new profitable avenues for operators. Furthermore, improvement in customer experience would also generate sizeable returns. The requirements laid out for 5G performance include enhancing the capacity by 1000 times, data rate by 100 times, a latency of less than a millisecond and 1000 times lower energy consumption. These requirements imply that significant challenges in the design of network architecture and air interface lie ahead.
In mobile communications, air interface is the crux of all the novelty undertaken in the development of a system. The progression in designing air interface schemes for 5G systems has taken considerable time over the years. Currently, filtered-OFDM (f-OFDM), which facilitates flexible waveform configuration, has been tipped as the 5G air interface scheme [27]. Under normal OFDM as in LTE-A, a specific numerology is used for the whole bandwidth, which results in limited spectral efficiency [28]. Under f-OFDM, the bandwidth is divided into several subbands for various types of services with appropriate numerology, resulting in enhanced utilization of spectrum ( Figure 1). the multiplexer departure process as an r-stage Coxian process. The authors of [25] present the performance evaluation of voice traffic in mobile radio networks with Coxian distributed channel holding time. The authors of [26] propose a method of performance analysis for a continuous time Markovian arrival process at the ATM Adaptation Layer 2 (AAL2), where the Common Part Sublayer (CPS) packets get multiplexed at the AAL2 buffer before being sent to transmission buffer. Both the buffers are modelled as queues.

Air Interface of 5G
The primary motivation behind the development of new standards for 5G mobile systems is innovations in support of IoT traffic. 5G applications for IoT are expected to increase market space and open new profitable avenues for operators. Furthermore, improvement in customer experience would also generate sizeable returns. The requirements laid out for 5G performance include enhancing the capacity by 1000 times, data rate by 100 times, a latency of less than a millisecond and 1000 times lower energy consumption. These requirements imply that significant challenges in the design of network architecture and air interface lie ahead.
In mobile communications, air interface is the crux of all the novelty undertaken in the development of a system. The progression in designing air interface schemes for 5G systems has taken considerable time over the years. Currently, filtered-OFDM (f-OFDM), which facilitates flexible waveform configuration, has been tipped as the 5G air interface scheme [27]. Under normal OFDM as in LTE-A, a specific numerology is used for the whole bandwidth, which results in limited spectral efficiency [28]. Under f-OFDM, the bandwidth is divided into several subbands for various types of services with appropriate numerology, resulting in enhanced utilization of spectrum ( Figure 1).

Relay Node Design and Scheduling Scheme
The RN has been introduced by 3GPP in Release 10 documentation [29]. RN is designed for extending the cell coverage area. The RN appears as a low-power BS to the mobile device. The RN is wirelessly linked to devices and the BS. The cost of deploying RN is less than deploying femto-and picocells, as no additional infrastructure is needed for connecting RN to the BS. RN improves coverage when placed at positions with weak channel conditions or dead spots in the coverage area.

Relay Node Design and Scheduling Scheme
The RN has been introduced by 3GPP in Release 10 documentation [29]. RN is designed for extending the cell coverage area. The RN appears as a low-power BS to the mobile device. The RN is wirelessly linked to devices and the BS. The cost of deploying RN is less than deploying femto-and picocells, as no additional infrastructure is needed for connecting RN to the BS. RN improves coverage when placed at positions with weak channel conditions or dead spots in the coverage area.
In this work, a wireless, fixed, inband, layer 3 type RN is proposed as the data aggregation and multiplexing node. This type of RN can be designed to aggregate and multiplex data packets from a large number of IoT devices. The proposed inband RN operates on the same radio frequency band that is specified to the BS by the network service provider. For accomplishing effective incorporation of IoT traffic in the 5G system, functionality of RN is modified slightly in this work.
Additionally, an RN Medium Access Control (MAC) packet scheduler [30] is used in this work to allocate frequency resources to IoT devices for RN access. The scheduler is capable of allocating resources to devices with various Quality-of-Service (QoS) traffic classes. In this work, the scheduling scheme used is blind equal throughput. This scheduler works in addition to the MAC layer scheduler placed at the BS for allocating resources to mobile devices and RNs connected to BS. It is also assumed in this work that all the IoT devices are having similar QoS requirements, therefore service differentiation is not performed. The BS scheduler in [31] is also employed in this work at the BS.
A multiplexing algorithm at RN is developed for aggregation of IoT traffic from different devices where packets are multiplexed before transmission to the BS. The multiplexing algorithm works collaboratively with the scheduler at RN to ensure that the size of multiplexed data is in accordance with the capacity of available radio resources ( Figure 2). In this work, a wireless, fixed, inband, layer 3 type RN is proposed as the data aggregation and multiplexing node. This type of RN can be designed to aggregate and multiplex data packets from a large number of IoT devices. The proposed inband RN operates on the same radio frequency band that is specified to the BS by the network service provider. For accomplishing effective incorporation of IoT traffic in the 5G system, functionality of RN is modified slightly in this work.
Additionally, an RN Medium Access Control (MAC) packet scheduler [30] is used in this work to allocate frequency resources to IoT devices for RN access. The scheduler is capable of allocating resources to devices with various Quality-of-Service (QoS) traffic classes. In this work, the scheduling scheme used is blind equal throughput. This scheduler works in addition to the MAC layer scheduler placed at the BS for allocating resources to mobile devices and RNs connected to BS. It is also assumed in this work that all the IoT devices are having similar QoS requirements, therefore service differentiation is not performed. The BS scheduler in [31] is also employed in this work at the BS.
A multiplexing algorithm at RN is developed for aggregation of IoT traffic from different devices where packets are multiplexed before transmission to the BS. The multiplexing algorithm works collaboratively with the scheduler at RN to ensure that the size of multiplexed data is in accordance with the capacity of available radio resources ( Figure 2). Thus, rather than demanding frequency resources from the BS for distinct IoT devices, resources are demanded by the RN for a bunch of IoT nodes. The IoT packets are aggregated and multiplexed at the RN. This multiplexing is achieved by, firstly, aggregating packets from several IoT devices and then multiplexing them into one large packet. On the BS side, a resource demand for the multiplexed large packet would be considered as one distinct radio resource demand. Therefore, sharing of radio resources among several devices would result in the enhancement of spectral efficiency. The capacity of the network would also increase, resulting in minimization of chances of network failure.

Data Aggregation and Scheme for Multiplexing
In this work, the design of a fixed RN is presented for 5G networks. IoT data traffic from devices located near RN is aggregated at the RN. In this paper, it is supposed that the RN antennas for the access link (interface with devices) and the backhaul link (interface with BS) are well separated and that self-interference can be totally avoided. For example, consider the case of deploying the access antenna inside some building while the backhaul antenna is deployed outdoors. In such cases, the Thus, rather than demanding frequency resources from the BS for distinct IoT devices, resources are demanded by the RN for a bunch of IoT nodes. The IoT packets are aggregated and multiplexed at the RN. This multiplexing is achieved by, firstly, aggregating packets from several IoT devices and then multiplexing them into one large packet. On the BS side, a resource demand for the multiplexed large packet would be considered as one distinct radio resource demand. Therefore, sharing of radio resources among several devices would result in the enhancement of spectral efficiency. The capacity of the network would also increase, resulting in minimization of chances of network failure.

Data Aggregation and Scheme for Multiplexing
In this work, the design of a fixed RN is presented for 5G networks. IoT data traffic from devices located near RN is aggregated at the RN. In this paper, it is supposed that the RN antennas for the access link (interface with devices) and the backhaul link (interface with BS) are well separated and that self-interference can be totally avoided. For example, consider the case of deploying the access antenna inside some building while the backhaul antenna is deployed outdoors. In such cases, the inband operation of RN can be ensured without any time division mechanism. At the MAC layer of the BS, a channel and service aware scheduler [31] is implemented for scheduling the RN as well as mobile user devices with regular data traffic. The scheduler consists of time and frequency domain uplink scheduling schemes.
An aggregation buffer is used at the RN for aggregating IoT data having a size that matches the size of instantaneously achievable Transport Block Size (TBS) offered by the air interface. The maximum buffer size denoted by n max corresponds to instantaneously available TBS, which depends on channel conditions. However, the multiplexed data forming a large packet can have a maximum size of n max − overhead (from lower layers). The overhead of LTE-A lower layers is assumed to be 352 bits [32]. The multiplexed large packet is transmitted to the BS via the backhaul link. This approach can bring significant improvement in the radio resource utilization, as discussed later. Conversely, this approach also can introduce problems due to constraints related to latency requirements associated with high priority data traffic, for example the delay sensitive emergency warnings. The small packets from IoT devices must wait until the buffer size reaches n max − overhead. In case of a highly loaded traffic scenarios, this issue may not occur, as the arrival rate is high and filling the buffer with packets would not take long. The issue arises in scenarios having low traffic load, where large interarrival times can cause longer waiting times greater delays in stuffing the buffer. As a result, the performance of delay sensitive IoT applications is compromised. This issue is dealt with by deploying a timer in the scheme. The timer is set to a certain time T max , which is the maximum waiting time of the first packet arriving since the empty buffer ( Figure 3). Thus, the buffer multiplexes the aggregated packets until a duration T max of the first arrival. Once waiting time of a packet reaches the timer expiry duration T max , the process of packet multiplexing is triggered, and the large multiplexed packet is sent to the BS. The original small packets from the IoT devices are then sent from the BS to the core network. inband operation of RN can be ensured without any time division mechanism. At the MAC layer of the BS, a channel and service aware scheduler [31] is implemented for scheduling the RN as well as mobile user devices with regular data traffic. The scheduler consists of time and frequency domain uplink scheduling schemes. An aggregation buffer is used at the RN for aggregating IoT data having a size that matches the size of instantaneously achievable Transport Block Size (TBS) offered by the air interface. The maximum buffer size denoted by corresponds to instantaneously available TBS, which depends on channel conditions. However, the multiplexed data forming a large packet can have a maximum size of − ℎ (from lower layers). The overhead of LTE-A lower layers is assumed to be 352 bits [32]. The multiplexed large packet is transmitted to the BS via the backhaul link. This approach can bring significant improvement in the radio resource utilization, as discussed later. Conversely, this approach also can introduce problems due to constraints related to latency requirements associated with high priority data traffic, for example the delay sensitive emergency warnings. The small packets from IoT devices must wait until the buffer size reaches − ℎ . In case of a highly loaded traffic scenarios, this issue may not occur, as the arrival rate is high and filling the buffer with packets would not take long. The issue arises in scenarios having low traffic load, where large interarrival times can cause longer waiting times greater delays in stuffing the buffer. As a result, the performance of delay sensitive IoT applications is compromised. This issue is dealt with by deploying a timer in the scheme. The timer is set to a certain time , which is the maximum waiting time of the first packet arriving since the empty buffer ( Figure 3). Thus, the buffer multiplexes the aggregated packets until a duration of the first arrival. Once waiting time of a packet reaches the timer expiry duration , the process of packet multiplexing is triggered, and the large multiplexed packet is sent to the BS. The original small packets from the IoT devices are then sent from the BS to the core network.

Analytical Model
The small IoT data packets from IoT devices arrive at the RN multiplexer. The first packet arrival at the multiplexer activates the timer. Due to the timer activation, the first packet would stay in the buffer for a duration not more than . During this time, arrival of additional packets would result in improvement of multiplexing gain. Once the packet waiting duration reaches the timer expiration threshold , the course of packet multiplexing is triggered. The large multiplexed packet is then transmitted to the BS.
is the maximum number of resource blocks available for RN to use in a single Transmission Time Interval (TTI). If some or all of the radio resource block are not needed in a TTI for RN traffic, the unallocated blocks are made available for regular mobile users

Analytical Model
The small IoT data packets from IoT devices arrive at the RN multiplexer. The first packet arrival at the multiplexer activates the timer. Due to the timer activation, the first packet would stay in the buffer for a duration not more than T max . During this time, arrival of additional packets would result in improvement of multiplexing gain. Once the packet waiting duration reaches the timer expiration threshold T max , the course of packet multiplexing is triggered. The large multiplexed packet is then transmitted to the BS. N max is the maximum number of resource blocks available for RN to use in a single Transmission Time Interval (TTI). If some or all of the N max radio resource block are not needed in a TTI for RN traffic, the unallocated blocks are made available for regular mobile users waiting for radio resources from the BS. Alternately, if the timer expiry time T max has not been reached, but the total size of small packets arrived at the multiplexing buffer reaches the size of n max − overhead; the large multiplexed packet is formed then and sent to BS before the timer expires. The size of n max depends on the available TBS, which in return is dependent on backhaul air interface radio channel conditions. At most, r packets can be aggregated into a single large packet and r is expressed as: where l describes the fixed size of the arriving packets in bits. Consider x as the actual number of packets in the RN buffer. As long as the waiting time does not reach T max or the buffer size does not reach n max − overhead, x remains in the interval 1 ≤ x < r. Packet arrival rate from IoT devices, λ is considered to be exponentially distributed. It can be shown that the total waiting time of x th packet is Erlangian [24]. Therefore, the probability of subsequent arrival of packet prior to multiplexing with x packets already in the buffer is denoted as a x and given as: The probability of initiating the multiplexing after x arrivals is b x where b x = 1 − a x . This multiplexer now emulates an r-stage Coxian process (Figure 4). The probability of an arrival with empty buffer, a 0 is 1 as the packet multiplexing starting without a single packet in the buffer is not possible. For triggering the multiplexing process, a single packet must arrive at the multiplexing buffer. Consequently, b 0 would always be 0. Upon arrival of the first packet from IoT device, the multiplexing system goes to stage 1 of the Coxian process with a total of r stages. Stage 1 illustrates that one packet has arrived in the buffer. During the first stage, two cases are possible. It is likely that the second packet also arrives into the multiplexing buffer ahead of timer expiry. The probability of this case is denoted as a 1 . But the converse case can also happen where no packet arrivals occur after the first one until timer expiry and the packet multiplexing is initiated with a single packet in the buffer. The probability of this case is expressed as b 1 . Now, if the second packet arrival occurs before the timer expires, the process is said to have entered stage 2. After entering stage 2, again two cases are possible. Either the third packet arrival occurs, which is considered to have a probability a 2 . Or the timer expires after the second arrival with a likelihood of b 2 . These probabilities a x and b x in each stage could happen only until x becomes equal to r, i.e., r packets arrive. As discussed earlier, r is the highest possible number of small packets that the buffer holds until starting the multiplexing process. Upon entering stage r, the sole probable event that can happen then is that the process of packet multiplexing begins immediately, and no more packet arrivals are awaited. At stage r, a r = 0 while b r = 1. the timer expires after the second arrival with a likelihood of . These probabilities and in each stage could happen only until becomes equal to , i.e., packets arrive. As discussed earlier, is the highest possible number of small packets that the buffer holds until starting the multiplexing process. Upon entering stage , the sole probable event that can happen then is that the process of packet multiplexing begins immediately, and no more packet arrivals are awaited. At stage , = 0 while = 1. The packet multiplexing process modelling is particularly useful in cases of data traffic where packet arrivals occur in bursts in between relatively silent periods. The generation of traffic at IoT devices is statistically independent. Thus, the traffic arriving at RN would be like bursts followed by periods of comparative silence. The goal of multiplexing is to combine the incoming packets at RN and make large blocks of data rather than sending isolated packets. This arrangement can reduce the The packet multiplexing process modelling is particularly useful in cases of data traffic where packet arrivals occur in bursts in between relatively silent periods. The generation of traffic at IoT devices is statistically independent. Thus, the traffic arriving at RN would be like bursts followed by periods of comparative silence. The goal of multiplexing is to combine the incoming packets at RN and make large blocks of data rather than sending isolated packets. This arrangement can reduce the number of blocks sent to BS considerably. The number of resource blocks needed for each stage can be determined by checking the TBS value for different resource block allocations and Modulation and Coding Scheme (MCS) combinations. Thus, the resource blocks required for different buffer sizes at the time of multiplexing is derived from the TBS table of 3GPP. The TBS values for M = 16 (as per 3GPP specifications [33]) and the buffer capacity in terms of data without overhead (352 bits) are given in Table 1 below. At stage 1, the buffer size in bits can reach l bits, where l is the fixed size of the arriving packets in bits. In this way, at stage x, the size of buffer would become x × l bits. The number of resource blocks needed for each stage x is denoted as N x,M , which relies on stage x and MCS M. The values of N x,M for M = 16 and l = 232 bits are given as in Table 2 below. To evaluate the performance of the proposed mechanism analytically, a discrete Random Variable (RV) X is defined, which denotes the number of packets fitted into a single multiplexed packet. The RV X can have possible values x = 1, 2, 3, . . . , r. Here, r can also be defined as the maximum number of IoT data packets inside a single multiplexed large packet. P[X = x] denotes the probability mass function of RV X. The likelihood of having more than x packets multiplexed into a large packet is a x . So a x is now defined in another manner. The probability a x represents the chance of having a value of X greater than x. This implies that the chance of getting more than x IoT data packets inside a multiplexed large packet is a x , which is given in the form of P[X > x] for 1 ≤ x < r as: This is illustrated with an example. Consider that a 1 = 0.6, which implies that 60% of the multiplexed large packets should be having a size greater than one small packet, which indicates that P[X > 1] = 0.6. Now the likelihood of X > 0 with x having value of 0 must always be 1, i.e., P[X > 0] = 1. This is due to the fact that multiplexing requires a minimum of one packet in the multiplexing buffer, which means that the value of X must be above 0 if multiplexing has to occur. Now by setting x = 0 in Equation (3), P[X > 0] turns out to be 1. In this fashion, the chance of X greater than or equal to r, i.e., P[X ≥ r] is always 0.
Since, X is a RV of discrete values, the likelihood P[X = x] is established as: This gives P[X = x] as:

Spectrum Utilization
Using N x,M from Table 2 and P[X = x] from Equation (8), the average resource blocks used in each multiplexing process or the average resource blocks consumed by a large packet is denoted by k, and determined as: The mean number of resource blocks consumed by a single multiplexed large packet, k is dissimilar from the parameter, mean number of blocks utilized per TTI. The mean number of resource blocks per TTI are considerably decreased by employing the multiplexing scheme along with the timer. However, k increases with usage of packet multiplexing scheme. The greater the value of k, the better the resource utilization. The multiplexing of packets happens if a duration T max has elapsed. If the packet aggregation scheme is not used at the RN and packets are sent to BS without multiplexing, then the arriving packets would be sent to BS immediately at the beginning of the next TTI, instead of waiting for the timer to expire. For this reason, the average number of resource blocks for each large packet would increase in case of multiplexing with the timer. As a result, the total number of multiplexed packets sent from RN to BS would be reduced considerably. The lesser the number of multiplexed packets, the higher the value of k and the lesser the overhead data.

Multiplexing Gain
Multiplexing gain is achieved by increasing the resource block utilization in every multiplexing process. For determining the multiplexing gain, the average resource block usage k without multiplexing (represented by k no_mux ) as well as with multiplexing (represented by k mux ) are determined in accordance with Equation (9). Multiplexing gain in percentage is given as:

Simulation Model
The analytical results of the derived model for resource block utilization are compared with simulation results for resource block usage. The simulation results in this paper are generated under simulation parameters given in Table 3. The evaluation underlines the case of relaying without multiplexing and the impact of relaying with multiplexing as well as timer mechanism. The two cases are compared to determine the multiplexing gain achieved both in analytical and simulation models. The OPNET simulation environment ( Figure 5) is used for determination of resource block usage and multiplexing gain. The OPNET model developed in this work emulates the impact of fast fading, slow fading, path loss, and noise on transmitted signal. It also models the interference from neighboring BSs. Channel models and interference models are developed according to [34]. The frequency reuse factor is 1. Therefore, devices transmitting uplink signals over certain resource blocks would cause interference in the transmission of devices of the neighboring BSs using the same resource blocks. The OPNET simulation environment ( Figure 5) is used for determination of resource block usage and multiplexing gain. The OPNET model developed in this work emulates the impact of fast fading, slow fading, path loss, and noise on transmitted signal. It also models the interference from neighboring BSs. Channel models and interference models are developed according to [34]. The frequency reuse factor is 1. Therefore, devices transmitting uplink signals over certain resource blocks would cause interference in the transmission of devices of the neighboring BSs using the same resource blocks.
where is path loss in dB and is the distance between receiver and transmitter in km. For slow fading, the log-normal model with standard deviation of 10 dB, mean 0 and correlation of 1 over distance is used. For fast fading, a time-dependent random function known as Jakes' model European Telecommunications Standards Institute (ETSI) based path loss models for different environments are used for research purposes. ETSI vehicular test environment model has been used in this work for modelling of path loss. L = 128.1 + 37.6 log 10 (R) (11) where L is path loss in dB and R is the distance between receiver and transmitter in km. For slow fading, the log-normal model with standard deviation of 10 dB, mean 0 and correlation of 1 over distance is used. For fast fading, a time-dependent random function known as Jakes' model is used. Thus, maps for coverage areas (cells) are generated using models in [34] and traces of these maps are fed into our OPNET simulation model. Noise and noise figure are modelled as per 3GPP specifications. Interference is modelled with in such a way that signals arriving from neighboring cells after reduction of power due to path loss, slow fading, and fast fading cause interference.
The RN OPNET model for relaying node is implemented to feature two radio interfaces, the Uu and the Un interfaces ( Figure 6). The corresponding protocol stacks are implemented as per 3GPP end-user protocols. The relaying scheduler is implemented in the MAC layer of RN Uu interface. The proposed packet aggregation and multiplexing scheme is implemented at the Packet Data Convergence Protocol (PDCP) layer of RN Uu interface. The GPRS (General Packet Radio Services) Tunneling Protocol (GTP) is also implemented where packet are tunneled at the RN and detunneling is performed at BS.  In the simulation model, the multiplexing transition probabilities are figured out by using a method of packet counting where is the likelihood of + 1 th packet arrival prior to the start of packet multiplexing such that packets are already present in the buffer and 1 ≤ < . In the simulation model, for determining the values of , is considered as the maximum number of small data packets in the buffer multiplexed to form a large packet. So could get a value in the range 1 ≤ < .
The analytical probability, is defined in (12). To find the a-posteriori probability for In the simulation model, the multiplexing transition probabilities are figured out by using a method of packet counting where a x is the likelihood of (x + 1) th packet arrival prior to the start of packet multiplexing such that x packets are already present in the buffer and 1 ≤ x < r. In the simulation model, for determining the values of a x , r is considered as the maximum number of small data packets in the buffer multiplexed to form a large packet. So x could get a value in the range 1 ≤ x < r.
The analytical probability, a x is defined in (12). To find the a-posteriori probability a x for simulation model, a record of the occurrences of a certain number of IoT packets packed inside a large packet is kept. For this purpose, a counter is deployed for recording arrivals at every stage of the aggregation buffer x. This means that a record of the number of instances where the Coxian process goes into a specific stage x is maintained. The counters designed for this purpose are denoted as c x for 1 ≤ x < r stages and defined in such a way that c x counts the number of times the large aggregated and multiplexed packet contained x packets or more. Thus c x is also the frequency instances that a particular stage x has been reached by the Coxian process.
Probability of a x is found with the help of value at counter c x . For example, to determine the probability a 1 , the value at counter c 2 is considered. This figures out the likelihood of arrival of second packet with respect to value at counter c 1 (i.e., the total number of arrivals). So, a x is calculated as:

Results and Analysis
In terms of k no_mux and k mux , comparison of analytical and simulation results for various traffic load scenarios are performed. The evaluation highlights the case of relaying without multiplexing and the impact of relaying with multiplexing plus timer mechanism. The performance of the multiplexing scheme is evaluated for scenarios with 100, 200, . . . , 900 up to 1000 IoT devices connected to the RN. These scenarios offer traffic load of less than 0.1. The reason for evaluating performance only in low load scenarios is that the expiry of timer rarely occurs in high load scenarios. The multiplexing process would almost always start before the timer expiry. The maximum load in the simulated scenarios is 0.07 for 1000 IoT devices in the RN coverage area.
Results are depicted using spider web charts. Such a chart is used to illustrate outcomes of multiple load scenarios in a single figure. Each chart axis represents one of 10 scenarios (i.e., 100, 200, . . . , 1000 users). Each axis signifies a scale of 1 to 5 ( Figure 7) for resource block utilization and 1 to 50 ( Figure 8) for multiplexing gain percentage. Figure 7 illustrates two aspects. First, the comparison of results for relaying with as well as without multiplexing depicts that the utilization of frequency resources is considerably improved. In the no multiplexing case, the average resource usage k is low. This shows that the total of multiplexed packets sent to BS from RN would be quite large. Therefore, the overhead data would grow significantly. The results for multiplexing with timer scenarios show that since packets are multiplexed into larger blocks, overheads needed for transmission blocks are significantly reduced. The second aspect is that the analytical results and simulation results agree with each other quite well. This conformity between analytical and simulation results indicate the validity of models. Simulation models are developed on the basis of artificial Random Number Generators (RNGs) for generation of random sequences of random values. Hence, simulation runs with various seeds for RNGs are performed to get results from varying random sequences. The 95% confidence intervals for simulation results are provided here in Table 4. only in low load scenarios is that the expiry of timer rarely occurs in high load scenarios. The multiplexing process would almost always start before the timer expiry. The maximum load in the simulated scenarios is 0.07 for 1000 IoT devices in the RN coverage area.
Results are depicted using spider web charts. Such a chart is used to illustrate outcomes of multiple load scenarios in a single figure. Each chart axis represents one of 10 scenarios (i.e., 100, 200, …, 1000 users). Each axis signifies a scale of 1 to 5 (Figure 7) for resource block utilization and 1 to 50 ( Figure 8) for multiplexing gain percentage.   Figure 7 illustrates two aspects. First, the comparison of results for relaying with as well as without multiplexing depicts that the utilization of frequency resources is considerably improved. In the no multiplexing case, the average resource usage is low. This shows that the total of multiplexed packets sent to BS from RN would be quite large. Therefore, the overhead data would grow significantly. The results for multiplexing with timer scenarios show that since packets are multiplexed into larger blocks, overheads needed for transmission blocks are significantly reduced. The second aspect is that the analytical results and simulation results agree with each other quite well. This conformity between analytical and simulation results indicate the validity of models. Simulation models are developed on the basis of artificial Random Number Generators (RNGs) for generation of random sequences of random values. Hence, simulation runs with various seeds for RNGs are performed to get results from varying random sequences. The 95% confidence intervals for simulation results are provided here in Table 4. The achieved multiplexing gain for all the traffic load scenarios is shown in Figure 8. Significant gains are achieved in all scenarios. Thus, it can be concluded that the framework provides better spectral utilization when IoT data traffic is sent over the network. It can also be noticed that the results for both simulation and analytical model agree with each other in Figure 8.   The achieved multiplexing gain G for all the traffic load scenarios is shown in Figure 8. Significant gains are achieved in all scenarios. Thus, it can be concluded that the framework provides better spectral utilization when IoT data traffic is sent over the network. It can also be noticed that the results for both simulation and analytical model agree with each other in Figure 8.
The M2M data traffic aggregation and multiplexing scheme introduces relaying delay to the uplink data packets. The evaluation of relaying delay is achieved by considering 3GPP MCS index 26, which implies relatively good channel conditions. The simulations are performed for two categories of environment settings. In one category, the simulations are performed without relaying and in the second category, simulations runs are carried out with RN for aggregation and multiplexing. The simulation parameters in addition to Table 3 parameters are provided here in Table 5. The simulations are performed for 15 scenarios, where the traffic load in each scenario varies from 1000 IoT devices up to 15,000 IoT devices. In the first scenario, the number of M2M devices is 1000. In each successive scenario, the IoT devices are increased by 1000. The load of data traffic, which is the ratio of arrival rate and service rate, is given in Table 6 for each scenario with various number of M2M devices. The packet delay performance of the system in both categories is compared in Figure 9. Error bars in the figure depict standard deviation. In low load scenarios, the delay performance without relaying is better than with relaying. The logic behind this behavior is that RN has to wait until arrival of n max − overhead bytes at RN, while no such waiting time is required if there is no relaying. However, once the traffic load increases, the environment without relaying fails to handle the load and large average delays (that do not fit to the scale used in Figure 9) are observed. In the scenario of 7000 IoT devices, this notion is clearly illustrated. However, in the environment with relaying, the traffic load of even up to 14,000 IoT devices is efficiently handled.

Conclusions
This paper provided a method for employing a low-cost 3GPP standardized layer-3, inband, RN for integration of IoT data traffic into 5G network. An RN design for aggregating packets and multiplexing IoT data before sending to BS was investigated in this work. An analytical model for data aggregation and multiplexing scheme was presented to determine resource block usage and multiplexing gain mathematically. The objective of the simulation model was justification of results achieved via analytical model using the arithmetic method. Results from both analytical and simulation models agreed with each other. Additionally, the simulation results for relaying delay performance were also illustrated. The analytical model proposed in this work provided a base for further development of analytical model for determination of transmission delays. The system performance evaluation with diverse traffic of various classes can also be derived in future.

Conclusions
This paper provided a method for employing a low-cost 3GPP standardized layer-3, inband, RN for integration of IoT data traffic into 5G network. An RN design for aggregating packets and multiplexing IoT data before sending to BS was investigated in this work. An analytical model for data aggregation and multiplexing scheme was presented to determine resource block usage and multiplexing gain mathematically. The objective of the simulation model was justification of results achieved via analytical model using the arithmetic method. Results from both analytical and simulation models agreed with each other. Additionally, the simulation results for relaying delay performance were also illustrated. The analytical model proposed in this work provided a base for further development of analytical model for determination of transmission delays. The system performance evaluation with diverse traffic of various classes can also be derived in future.


Intervention Model of Sinabung Eruption to Occupancy of the Hotel in Brastagi, Indonesia

This paper presents an analysis of intervention of Sinabung eruption to the rate of hotel reservation in Brastagi, North Sumatera, Indonesia. The weekly data of this study was collected from four hotels in Brastagi for the period of January 2012 – May 2015. The method that be used is an analysis of the intervention with ARIMA (Autoregressive Integrated Moving Average) process by Box – Jenkins. To see the stationary of data is used the plot of ACF (Autocorrelation Function). For identification of model is used the plot of PACF (Partial Autocorrelation Function), and to choose the best of model is used AIC (Akaike Information Criterion). Meanwhile to see the suitability of data to model is used Q-Statistics Ljung-Box. The result shows that decline the rate of hotel reservation only until the 30th period after the intervention. After that, the rate of hotel reservation will be as normal as before the eruption. The results of this study are limited, valid only on the subject of research, not generally applicable, to be generalizable needed further research.

INTRODUCTION:
The eruption of Mount Sinabung in Karo district, province of North Sumatera, Indonesia on September 15, 2013 had a negative impact on the tourism sector in Brastagi city which is approximately 11 km from Mount Sinabung, especially in term of hotel reservation in the surrounding areas impacting the city's income. Intervention in time series is an event or events that may affect the data pattern of the time series. Intervention model in time series analysis was introduced by Box and Tiao in 1975 that discussed the effect of the enactment of legislation engine design to the level of pollution oxidant in Los Angles (Box & Tiao, 1975). Others the paper related to this research are (Montgomery & Weatherby, 1980) investigated the influence of the Arab oil embargo against the level of electricity consumption in the United State. (Suhartono & Huriroh, 2003) investigated the influence of New York WTC bombing against fluctuations in the price of the world. (Chung, Ip, & Chan, 2009) using ARIMA intervention model in research, investigated the effects of the financial crisis on manufacturing industry in China, while (Jarrett & Kyper, 2011) also uses the ARIMA intervention model to predict and analyze the Chinese stock prices. This paper focuses on the determination of the intervention model of Sinabung eruption to the rate of hotel reservation in Brastagi city.

MATERIAL AND METHODS:
The data used in this study is the rate of hotel reservation in Brastagi based on the number of rooms booked by visitors. Data collected weekly from four hotels in Brastagi, two hotels located alongside Sumatera highway, namely Grand Mutiara hotel and The Hill Sibolangit, while two more hotels are located alongside the road to Gundaling, namely Sibayak International Hotel, and Brastagi Cottage. All four hotels are seen to represent all hotels in Brastagi. This study does not consider the peak time, because the data is weekly and the peak time occurs on weekens. Sampling strategy using cluster sampling for determining the location, then the selection of hotels is done randomly. The number of time series data the hotel reservation rate there are 178 of data, starting on the first week in January 2012 as the 1 st data until the fourth week of May 2015 as the 178 th data. Based on the data used, the first Sinabung eruption at second week of September 2013 viewed as the intervention on the 89 th data. The method used is the analysis of intervention with ARIMA process.

RESULTS:
The data ploting when eruption ocucured at t = 89, as Figure 1 below:

Figure 1: Plot of time series data and Intervention Analisys before Intervention:
Based on the Sinabung eruption on the second week of September 2013 viewed as the intervention on the 89 th data, therefore there are 88 datas before the intervention on the 1 st data through to 88 th data and are denoted as x1. The data plotong, ACF, and PACF consecutive for x1 presented in Figure 2, 3, and 4 below: ACF plot in Figure 3 in the form damped sine wave, this indicates that the time series data for x1 is stationary. ACF plot in Figure 3 is disconnected after first lag and PACF plot in Figure 4 is also disconnected after first lag, this indicates that the model is appropriate for x1 series is ARIMA (1,0,0), ARIMA (0,0,1) and ARIMA (1,0,1). Furthermore, the calculation with the software "R", obtained the value of parameters and the value of AIC as Table 1 below: To conclude whether the model of equation (1)  Testing with Q-statistic Ljung-Box: Chi-squared = 0.0069, df = 1, p-value = 0.9337. Value statistics count < value of Statistics Tables, so it can be concluded that the model is fit data.

Analisys of Intervention:
Identification of intervention response which has been done by observing the data pattern before, during, and after the intervention. Observations were made on the data pattern around the 89 th data, ie from the 70 th data until the 108 th data, presented in Figure 6  In Figure 6 can be seen that the effect of interventions at the 2 nd week of September 2013 (89 th data). Although not decline dramatically, however, after the eruption seen a decline until some lag time in the future, so that the function of intervention for this case is: model is appropriate for series x2 is ARIMA (1,0,0), ARIMA (0,0,1), ARIMA (1,0,1), ARIMA (0,0,2), ARIMA (2,0,0), ARMA (1,0,2), ARMA (2,0,1) and ARMA (2.0, 2). Furthermore, the calculation with the software "R", obtained the value of parameters and the value of AIC as Table 2 below: (3) To conclude whether the model of equation (3) fits to the data, used the autocorrelation test for residuals by Q-Statistic Ljung-Box, to test whether the residuals   t a is a sequence of random variables with mean is zero and constant variance, as Figure 10 below: Testing with Q-statistic Ljung-Box: chi-squared = 0.0283, df = 1, p-value = 0.8665. Value of statistics count < value of Statistics

CONCLUSION:
Plot of the time series data in figure 11 shows that after the Sinabung eruption on 15 September or the 3 rd week of September 2013 as the 89 th data, the data pattern changes that inflict decrease of hotel reservation until the 119 th data. After the 119 th data, the pattern back into its previous form before the intervention, this indicates that the imfact of intervention on the hotel reservation only until the 119 th data or for 30 weeks after the eruption, after that, the rate of hotel reservation will be as normal as before the eruption. Figure 11: Limit the impact of intervention


In Search for the Super Element: Algorithms to Generate Higher Order Elements

This work was motivated by the need that may arise during the creation of finite element programs, for higher order elements. The problem of selecting or creating shape functions that satisfy the required need of the problem may be one major problem that stand in the way of the element creation. In this work we present an attempt that points in the direction of creating the element matrices for higher order elements using simple, Lagrange, and modified Lagrange polynomials. The results obtained for the test cases indicate the possibilities and limitations on those attempts. It is concluded that the modifies Lagrange polynomials have a great potential for use in elements with high number of nodes and/or high number of DOF per node. Nevertheless, they impose a high computational cost on the element matrices generation.

INTRODUCTION
In the past decades, the finite element method has grown into a well developed numerical technique that can be used to provide accurate and relatively quick solution to boundary value problems that describe several physical problems. The procedure of the finite element model starts by dividing the domain into elements that are connected through nodes. Each node has degrees of freedom (DOF) that, usually, present the value and derivatives of the physical function of interest. Then the elements are defined by the interpolation functions that describe the change of the function within their boundaries. Following that, the whole domain description is created by assembling the elements and applying the boundary conditions. Finally, the problem is solved to obtain the values of the DOF of all the nodes. 1,2 The basic assumption is that as the number of elements increase in the domain, the solution approaches the exact one. The power of the finite element method lies in that it may use simple functions, usually polynomials, to describe the change of the target function inside an element, and then use that to describe the change in the whole domain by assembling the different elements that cover the domain.
Another approach to reach higher accuracy is to increase the order of each of the elements used in the model which reduces the need for higher number of elements (p-version finite element). Higher order elements may be generated by using more nodes in the element, using more complex functions, or by hierarchical elements. 3,4,5,6 Most of the work performed for the solution of the finite element problem is algorithmic, that enabled the excessive use computers and the production of computer packages that can manipulate several physical problems with complex domain geometries. Nevertheless, in the core of all finite element codes, lies the element matrices. The element matrices are generated through the knowledge of the mathematical model that describes the physical problem and the selection of the interpolation functions, shape functions. Equipped with those two pieces of information, you may generate the element matrices for any physical problem. However, the generation of a general element is still far from algorithmic. 7 A special problem faces the generation of the finite element model is the selection of the interpolation function that satisfies the prescribed degrees of freedom at the nodes 5 When the DOF involve derivatives of the function, the Lagrange polynomials should be selected to satisfy the condition that they have the slope of one at the associated node and zero at the others as well. However, it is quite difficult to find the Lagrange polynomials that may be used in a general problem that is why they may be generated by using classical interpolation techniques.
The classical interpolation techniques require the evaluation of the polynomial coefficients by forcing the polynomial to satisfy the boundary conditions, thus, requiring the solution of a set algebraic equation, hence, the inversion of a matrix. The matrix inversion process is a straight forward numerical technique; however, as the number of variables increase, the matrix eventually becomes singular because of the round-off errors. 9, 10 Another challenge that faces the creators of finite element models is the performance of integration. The element matrices are generated by integrating the shape functions or their derivatives over the element. In many simple cases, the integration is readily available by hand or using symbolic manipulators. However, in more practical cases, numerical integration is required which, in turn, introduces errors to the resulting element models.
In the past years, several research articles were published, mostly by mathematics oriented researchers, about the topic of generating higher order polynomials, shape functions, for applications in the finite element models 5, 6, 8, 11, 12. Such research, among many others not cited here, accomplished a great task in the direction of generalizing the automation of finite element model generation.
Following that, an excellent attempt for the automation was presented in 7 using the symbolic manipulator Mathematica® and an accompanying package AceGen®.
None of the research reviewed by the author presented a straight-forward methodology that may be used by the engineers who need to create finite element models on numerical manipulators, thus, in the following work, we will be attempting to describe a generalize procedure for the generation of the finite element matrices with an eye on elements with large numbers of nodes and degrees of freedom. We will start by describing a classical method of generating the matrices and creating the shape functions, then we will move towards describing attempts to generate a general method for performing exact integration for generating element matrices. Following that, we will describe the attempts using Lagrange polynomials to avoid matrix inversion, finally, modifies Lagrange polynomials will be created to enable to the creation of different finite element problems.
The main purpose of this work is providing a methodology that may by applied using numerical manipulators, especially open source, to generate the element matrices that may be used by researchers without the need for commercial packages or commercial symbolic manipulators.

CLASSICAL APPROACH 2.1. Interpolation Polynomials -Shape Functions
The finite element model starts by assuming the solution of the unknown function in terms of an interpolation polynomial that satisfies the values of the degrees of freedom at the nodes of the element. In most text, these functions are denoted the letter N such that: Where fi are the values of the function at node i and Ni(x) are polynomials that have the value of one at the corresponding node and zero at every other node. The above relation may be written as: The first problem that faces the researcher who wants to create a new finite element model is to find and select the polynomials that may be used for the problem. However, that is not a great problem for low order elements since the polynomials are readily available in many texts. Nevertheless, coding them may involve mistyping which will require debugging of the code. On the other hand, it becomes a bit inconvenient to change the number of nodes per element to test the effect of higher, or lower, order elements. Thus we may use another way of defining the polynomials that can reduce both problems, that is by using regular simple polynomials such that: Where ai are constants, generalized coordinates that should be determined for the function to satisfy the given values. To do that we will need to set a number of equations and solve them to find those values, thus: Which may be written as: {δ } It can be readily proven that: Thus using this procedure we were able to obtain the polynomials for any number of degrees of freedom without having to look them up in literature or enduring the problems of mistyping them into the code. However, nothing comes for free; using this method will involve matrix inversion which will become ill-conditioned as the number of degrees of freedom increase.
If the degrees of freedom of the element involve the derivatives of the function, we may adjust the above procedure. Given the value of the function fi and the fist derivative f'i, for example, we may write: Then the [T] matrix may be constructed such that: In this case, the number of degrees of freedom per element will be equal to twice the number of nodes and the function Ni(x) will appear in pairs. The first function of the pair will have a value of one at its corresponding node while the slope at that node is zero, and has a value and slope of zero at every other node. Meanwhile, the second function of the pair will have a slope of one and a value of zero at the corresponding node while its value and slope are equal to zero at every other node.
Note that, in both cases above, the order of the polynomial terms in the H(x) matrix does not affect the resulting shape functions because they are going to be rearranged using the T matrix.

Element Matrices
Each finite element problem will generate a set of matrices for each element. These matrices are derived from the physical model that is usually presented in the form of a differential equation. The matrices are constructed by integrating some derivative of the function that usually appears in the form: Where Q(x) is some function that describes the physical properties of the problem, D m (.) is a differential operator that differentiates the function f(x) m times. For example, in mechanics of material, the finite element matrix for a bar in static problem may be obtained using: Where E(x) and A(x) are the modulus of elasticity and the cross-section area respectively. For simple problems, the above integral may be performed by hand or using symbolic manipulators. However, it is more convenient to use numerical integration in order to add flexibility to the programming. Nevertheless, some cases will not allow for hand or symbolic manipulators' integration, thus it is common to find numerical integration routines associated with any finite element program. One of the popular techniques used for numerical integration is the Gauss-Legendre quadrature which produces almost exact integration for polynomials that may reach up to 24 th order. 12

2-D Problems
In 2-D problems, the interpolation polynomials should include combinations of the x and y coordinates. To identify the terms, the common practice is to select them from what is known as the Pascal triangle (See Figure 1). The number of terms selected should be equal to the number of degrees of freedom and it is recommended that the terms should create a symmetric shape in that triangle.

Figure 1. Pascal triangle
The Pascal triangle is quite handy in many problems, however, if you strict your work to using the full polynomials in both directions, then the generation of the terms becomes quite straight forward by multiplying both x and y-polynomials. For example, for a quadrilateral element with a single degree of freedom per node, we need a linear polynomial in each of the directions, hence, we may write: →H ( x, y)= ⌊1 x y xy ⌋ Then we proceed with the H(x,y) vector just as in the 1-D cases: Note, again, that the order of the polynomial terms in the H vector will not affect the resulting shape functions similar to what we had in the 1-D problems. At this point, we still have both problems when generating the higher order elements, namely, the matrix inversion and the numerical integration. In the following section, we will attempt generating the matrices using exact integration to avoid the numerical integration problem.

EXACT INTEGRATION
To avoid the errors introduced by the numerical integration, we may resolve to exact integration using symbolic manipulators (wxMaxima® , Mathematica® , Mable® , etc …) but then we are back to the problem of rewriting the element matrices into our code, with all the problems that may include, or write the whole finite element model using the symbolic package which is normally extremely slow when it comes to numerical manipulations compare to numerical coding packages (Octave® , Matlab® , etc ...) or programming languages (Fortran, C++, etc …). Thus, we resolve to write an algorithm to perform the exact integration on the numerical package for specific elements with unknown number of degrees of freedom. This is where the H(x) row-matrix becomes very handy. Since the row-matrix is composed of simple polynomial terms, their differentiation, and later integration, is readily evaluated.
Let's examine the case of a bar element with n nodes (and n degrees of freedom). For the bar element, the element matrix, stiffness matrix, is evaluated using the integration: For the sake of the illustration, we will assume that the modulus of elasticity and the cross-section area are constants, thus, we may write: But, the transformation matrix is also constant, thus, it may be dragged out of the integration leaving us with the derivative of the H(x) vector.
In this case, the T matrix may be given as: If we examine the terms of the H(x) vector, we may write: Thus the first derivative may be written as: International Journal of Trend in Scientific Research and Development (IJTSRD) ISSN: 2456-6470 Page: 512 Hence, we may write the matrix: From that, we may get the integration: Similarly, for any 1-D problem involving any derivative, we may readily obtain the general term for the G matrix and perform the integration to obtain the G* matrix.
Let's examine the case of a beam element with n nodes (and 2n degrees of freedom). For the beam element, the element matrix, stiffness matrix, is evaluated using the integration: For the sake of the illustration, we will assume that the modulus of elasticity and the cross-section second moments of area are constants, thus, we may write: But, the transformation matrix is also constant, thus, it may be dragged out of the integration leaving us with the derivative of the H(x) vector.
In this case, the T matrix may be given as: If we examine the terms of the H(x) vector, we may write: Thus the second derivative may be written as: From that, we may get the integration:

, n
For 2-D problems, the problem becomes a little more complex. The general term of the H(x,y) vector will depend on the order at which we select the terms from the Pascal's triangle. For illustration, let's assume that we are going to generate a rectangular element for plate bending that ensures continuity of slope (C 1 element).
In this case, the number of nodes in the element may be selected to be n 2 where n is the number of nodes per side of the element (n=2,3,4, …) in this case the number of degrees of freedom per element will be 4n 2 , and the x and y polynomials, each, will be of 2n-1 order. We may write:

,2n
And the derivatives of the H vector, may be written as: The element stiffness matrix for the static loading problem of thin plates may be written as: Where: And: Whose general term may be written as: [Q] is the plate stiffness which may be written, for the sake of compactness, in the form: The matrix [G] may be written as: where the general term may evaluated by: Similarly: From which, we may get the integral of the general term of the [G] matrix in the form: a(α 1 +α 4 )+b(α 2 +α 3 )+4Cα 5 In both the 1-D and 2-D problems illustrated above, we were able to create the element matrices exactly using the trick of integrating the H vector and its derivatives. However, the limiting condition on generating high order elements will always be the inversion of the T matrix which will become ill-conditioned with higher order elements. For example, the plate bending element described above worked only until the number of nodes became 16 (64 DOF) (7 th order polynomial in each of the x and y directions) using double precision numbers on the Octave numerical manipulator. Hence, a search for another technique was required.

LAGRANGE POLYNOMIALS 4.1. 1-D problems
The Lagrange interpolation polynomial for an (n-1) th order polynomial, which may be used when n degrees of freedom are give, may be written in the form: Thus, we may obtain the first derivative in the form: We may also obtain further derivatives for use in the element equation, however, let's remember that the function interpolation is based on the values of the function, thus, if we need to have any of the derivatives of the function as a degree of freedom, we cannot use this interpolation method. Nevertheless, we have obtained the shape functions without having to go through the matrix inversion problem described in the previous section. On the other hand, performing the exact integration will not be readily available for the general term since it now involves multiplication of several linear terms. Luckily, as mentioned earlier, the Gauss quadrature integration provides high degree of accuracy for a relatively high order polynomial integration.
Note that he order of the Lagrange polynomials used above has to coincide with the order of the node numbering unlike the cases we described in sections 2 and 3 above.

2-D Problems
With the same advantages and disadvantages of the 1-D problem, the derivation of the shape functions for a 2-D problem comes straight forward by multiplying the polynomials in the x and y directions.
In a vector, we usually use a single index to point to the different element, thus, we need to create a relation between the vector index and the i and j counters of the nodes counters in a 2-D element. A simple relation for counting the nodes in the x-direction first would be in the form: Thus we may write the general shape function in the form: For which we may evaluate the first derivatives to be:

MODIFYING THE INTERPOLATION FUNCTION
Borrowing from the Spline interpolation techniques, we may work around with the Lagrange polynomials to obtain polynomials that may include derivatives of the function. In the following, we will work with functions that are required to satisfy the value of the function and its first derivative at each node of the element. We will start be deriving the cubic polynomial that will fit 2 values and two slopes, then, we will describe the process for quintic polynomial, followed by the general approach for obtaining the (2n-1) th polynomial that may be used for n-node elements.

The Polynomials
Using the similarities we obtained in the 3 rd and 5 th order polynomial derivation, we will present a generalization in the following section. We may write the general relation for a function of the (2n-1) th order that uses n pairs of polynomials as: Where: This function needs to satisfy the conditions: For the first condition, we get: To proceed with the second condition, we get the slopes as: Substituting: Which gives: Thus we may write the pair of polynomials as: Rearranging, we obtain: From which we may write the shape functions as: In 2-D problems, the interpolation functions will be a straight forward generalization of the above procedure giving the function in the form: n= n x × n y Where: 1,2,... ,n x j= 1,2,... ,n y m= 1,2,. ..,n This function needs to satisfy the conditions: Which will result in a quartet of shape functions in the form:

The Derivatives
In the problems that make use of the shape functions described in the previous section, the finite element model usually requires the first and second derivatives of the shape function to perform the required calculations. In this section, we will present those derivatives for the 1-D problem without elaborations. Table 1 presents the change of the values of the minimum Eigenvalue of the stiffness matrix generated for a bar element. If we ignore the numbers that have an absolute value less than 10 -10 considering them as a numerical zero, we may still find that increasing the number of nodes per element start to introduce negative Eigenvalues until we reach the failure criterion set in this work as 10 -4 . Note that, unpredictably, the exact integration method failed before the numerical integration at 15 nodes per element. On the other hand, the Lagrange polynomial method, which does not use the transformation matrix, was able to produce satisfactory results up to 23 nodes per element which approaches the limit of the 12 point numerical integration technique used. Table 2 presents the results obtained for a beam element. Again we find that the exact integration and numerical integration elements, both, failed around the 18 and 16 DOF respectively which coincides with the number of DOF at which failure occurred in the bar element. Meanwhile, the Modified Lagrange method was able to continue all the way up to 24 DOF which reaches the limit imposed by the numerical integrator.    Table 3 presents the results for a plate bending problem. The same patterns observed in the bar and beam elements can be observed in the plate problem for the modified Lagrange case.

CONCLUSIONS
Several problems were examined in this work to demonstrate the limitations imposed on the creation of super elements with high number of nodes and degrees of freedom. The element matrices were generated using simple, Lagrange, and modified Lagrange polynomials and the integration was performed numerically using Gauss-Lagrange quadrature or exact when possible. Problems involving 1-D elements with 1 or 2 DOF per node and 2-D problems with 4 DOF per node were examined by increasing the number of nodes per element until the lowest (negative) Eigen value of the stiffness matrix had an absolute value greater than 10 -4 . From the results obtained above we may conclude the following: ➢ Automated generation of higher order isoperimetric C 0 and C 1 continuous element model was enabled using an open source numerical manipulator (Octave®) using the different methods presented in this work. ➢ Generating higher order elements is computationally expensive compared to lower ones. Thus, the h-version is more efficient for the same accuracy than the p-version.
➢ Using numerical integration to evaluate the element matrices is efficient and accurate enough for most practical purposes ➢ Using Lagrange and modified Lagrange polynomial increases the stability of the results for higher order element problems when they become a necessity to use. In such cases, the limitation on the results' accuracy will be imposed by the accuracy of the numerical integration scheme.
Further work with the modified Lagrange polynomials, proposed in this work, is needed to enable the inclusion of higher order derivatives in the interpolation polynomials. Such research will help in the generation of elements with the needed derivatives without falling into the trap of the ill-conditioned transformation matrix.