• Author / Uploaded
• claudiunicola

• Categories
• Análisis de conglomerados
• Json
• Áreas de informática
• Informática
• Tecnología

Citation preview

Applications for the

MATLAB Practical Engineer Edited by Kelly Bennett

MATLAB Applications for the Practical Engineer Edited by Kelly Bennett

D3pZ4i & bhgvld, Dennixxx & rosea (for softarchive)

Stole src from http://avaxho.me/blogs/exLib/ Published by AvE4EvA Copyright © 2014 All chapters are Open Access distributed under the Creative Commons Attribution 3.0 license, which allows users to download, copy and build upon published articles even for commercial purposes, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications. After this work has been published, authors have the right to republish it, in whole or part, in any publication of which they are the author, and to make other personal use of the work. Any republication, referencing or personal use of the work must explicitly identify the original source. As for readers, this license allows users to download, copy and build upon published chapters even for commercial purposes, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications. Notice Statements and opinions expressed in the chapters are these of the individual contributors and not necessarily those of the editors or publisher. No responsibility is accepted for the accuracy of information contained in the published chapters. The publisher assumes no responsibility for any damage or injury to persons or property arising out of the use of any materials, instructions, methods or ideas contained in the book. Publishing Process Manager Technical Editor AvE4EvA MuViMix Records Cover Designer

Published September 08, 2014 ISBN-10 953511719X ISBN-13 978-9535117193

Contents Preface

Chapter 1 Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon Web Services (AWS) by Kelly Bennett and James Robertson Chapter 2 A MATLAB-based Microscope by Milton P. Macedo Chapter 3 Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application by R. Silva-Ortigoza, C. Mÿrquez-Sÿnchez, F. Carrizosa-Corral, V. M. Hernÿndez-Guzmÿn, J. R. Garcÿa-Sÿnchez, H. Taud, M. Marciano- Melchor and J. A. ÿlvarez-Cedillo Chapter 4 Dual Heuristic Neural Programming Controller for Synchronous Generator by Mato Miskovic, Ivan Miskovic and Marija Mirosevic Chapter 5 Design of Fractionation Columns by Hassan Al-Haj Ibrahim Chapter 6 Remotely Train Control with the Aid of PIC32 Microcontroller by Mostefa Ghassoul Chapter 7 Integration of MATLAB and ANSYS for Advanced Analysis of Vehicle Structures by A. Gauchÿa, B.L. Boada, M.J.L. Boada and V. Dÿaz Chapter 8 Modelling and Analysis of Higher Phase Order (HPO) Squirrel Cage Induction Machine by A.A. Jimoh, E.K. Appiah and A.S.O. Ogunjuyigbe Chapter 9 Atmospheric Propagation Model for Satellite Communications by Ali Mohammed Al-Saegh, A. Sali, J. S. Mandeep, Alyani Ismail, Abdulmajeed H.J. Al-Jumaily and Chandima Gomes Chapter 10 Stateflow® Aided Modelling and Simulation of Friction in a Planar Piezoelectric Actuator by G. M’boungui, A.A. Jimoh, B. Semail and F. Giraud Chapter 11 Modelling and Simulation Based Matlab/Simulink of a Strap-Down Inertial Navigation System’ Errors due to the Inertial Sensors by Teodor Lucian Grigorie and Ruxandra Mihaela Botez

VI

Contents

Chapter 12 Tool of the Complete Optimal Control for Variable Speed Electrical Drives by Marian Gaiceanu Chapter 13 Modeling of Control Systems by Roger Chiu, Francisco J. Casillas, Didier Lÿpez-Mancilla, Francisco G. Peÿa-Lecona, Miguel Mora-Gonzÿlez and Jesÿs Muÿoz Maciel Chapter 14 Advanced Decimator Modeling with a HDL Conversion in Mind by Drago Strle and Duÿan Raic Chapter 15 Code Generation From MATLAB – Simulink Models by Tiago da Silva Almeida, Ian Andrew Grout and Alexandre Cÿsar Rodrigues da Silva Chapter 16 DC/DC Boost-Boost power converter as a didactic system: Experimental results with Matlab/Simulink via current and voltage probes by R. Silva-Ortigoza, M. Antonio-Cruz, M. Marcelino-Aranda, G. Saldaÿa-Gonzÿlez, C. A. Merlo-Zapata and P. Pÿrez-Romero Chapter 17 Using Matlab and Simulink for High – Level Modeling in Biosystems by Cristina-Maria Dabu Chapter 18 Eigenvalue Analysis in Mode Domain Considering a Three-Phase System with two Ground Wires by R. C. Monzani, A. J. Prado, L.S. Lessa and L. F. Bovolato Chapter 19 Analysis of Balancing of Unbalanced Rotors and Long Shafts using GUI MATLAB by Viliam Fedÿk, Pavel Zÿskalickyÿ and Zoltÿn Gelvanic Chapter 20 Analysis of Robotic System Motion in SimMechanics and MATLAB GUI Environment by Viliam Fedÿk, Frantiÿek Durovskyÿ and Rÿbert ÿveges Chapter 21 Vibration Analysis of Laminated Composite Variable Thickness Plate Using Finite Strip Transition Matrix Technique and MATLAB Verifications by Wael A. Al-Tabey Chapter 22 Image Processing with MATLAB and GPU by Antonios Georgantzoglou, Joakim da Silva and Rajesh Jena

Preface

Being an automotive enthusiast and hobbyist, there is a saying around the race track that speed isnt everything, its just most everything. Similar sentiments can be said about MATLAB. MATLAB isnt everything to a practicing engineer, its just most everything. MATLAB is an indispensable asset for scientists, researchers, and engineers. The richness of the MATLAB computational environment combined with an integrated development environment (IDE) and straightforward interface, toolkits, and simulation and modeling capabilities, creates a research and development tool that has no equal. From quick code prototyping to full blown deployable applications, MATLAB stands as a de facto development language and environment serving the technical needs of a wide range of users. As a collection of diverse applications, each book chapter presents a novel application and use of MATLAB for a specific result.

Chapter 1

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon Web Services (AWS) Kelly Bennett and James Robertson Additional information is available at the end of the chapter http://dx.doi.org/10.5772/58895

1. Introduction The popularity of social networking has allowed access to staggering amounts of unique data, which has created new possibilities for data analysis and exploitation. Such data has proven useful in marketing, decision making, destabilizing terrorist networks, behavior evolution, and determining future social trends [1]. Increased usage of social networking sites has also been observed during events related to natural disasters; significant political, sporting, and social events; and other crises. Twitter users provide status updates through tweets. Since tweets are typically short and always less than 140 characters, they may need to undergo additional analysis to provide contextual clues. Applying traditional natural language processing algorithms on such data is challenging. Before making use of this data, it must first be extracted, processed and analyzed appropriately using certain algorithms and theories. According to a 2011 study from the International Data Corporation (IDC)–a marketing firm specializing in information technology and other consumer technologies–unstructured data, that is, data that does not have a pre-defined data model or is not organized in a pre-defined manner, is growing at a faster rate than structured data. Within the next decade, unstructured data will account for 90 percent of all data created. A large driving factor in the increase in unstructured data is social networking data, such as tweets. It is estimated that more than 80 percent of all potentially useful data is unstructured [2].

4

MATLAB Applications for the Practical Engineer

The success of businesses in the coming decade will likely rely on their ability to successfully analyze data from social networks. In addition to the increase in social networking, cloud computing technologies continue to increase in popularity and provide end users the ability to quickly spin up powerful servers with analytical capabilities, making analysis of large amounts of data more affordable and practical than in previous decades. Main topics addressed in this chapter include: • Explore viable approaches for extracting social networking data sets using tools compatible with MATLAB1 and cloud technologies and infrastructures. • Use existing data mining and statistical tools within MATLAB to conduct analysis on social networking site data. • Discuss potential cost savings and document approaches for implementing social network‐ ing site data analysis via the cloud through vendors such as Amazon Web Services (AWS). • Provide a tutorial on the use of MATLAB tools to analyze unstructured data with an emphasis on social networking data. • Provide an overview and example on the use of MATLAB in a commercial cloud environ‐ ment, such as AWS.

2. Previous work Muhammad Mahbubur Rahman, author of “Mining Social Data to Extract Intellectual Knowledge,” was able to extract data from Facebook using various data mining techniques. He first determined the most frequent terms that were being used, such as birthday, about me, gender, and music, and then created a database using these terms and the appropriate data types. [1] From the “about me” information, he was able to use an Euclidian distance formula to put each user into a certain “class”, for example, aggressive, dishonest, romantic, eager to learn, or lazy. Grouping classes by age and gender creates a method of comparing various attributes by age and gender. This type of research demonstrates the ability to use data mining techniques to provide intellectual knowledge that can be used to possibility predict human behavior, provide insight into human decision making, or provide a method of determining anomalous events. Data mining has proven useful in Twitter as well, with much research performed into detecting trends and bursty keywords. One has to pay attention to some abnormalities with Twitter, such as many people use smart phones to tweet and it is common for people to make mistakes while typing on smart phones. One would have to account for misspellings and merge them 1 Registered trademark of Mathworks Corporation.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

into one variable so they are not separated and lost. In addition, some phrases are commonly abbreviated, such as “NYC” for “New York City,” or reduced, such as “Vegas” for “Las Vegas.” A study, using algorithms to account for these abnormalities, was performed that determined trendy and bursty keywords on Twitter with 92% accuracy when compared with Google trends [3]. Karandikar [4] described an approach to determine the most suitable topic model to cluster tweets by analyzing the effect of change in topic model parameters such as training data size and type and the number of topics on its clustering performance. The model was able to cluster tweets, using this approach, with an accuracy of greater than 64 percent for the two specific events. Perera [5] developed a software architecture using a Twitter application program interface (API) to collect tweets sent to specific users. The extracted data were processed to characterize the inter-arrival times between tweets and the number of re-tweets. Analysis revealed that the arrival process of new tweets to a user can be modeled as a Poisson process while the number of retweets follows a geometric distribution. Turner and Malleson [6] presented an exploratory geographical analysis of a sample of Twitter post data from June 2011 to March 2012 in the city of Leeds, England. Geographical cluster detection methods were used in combination with text mining techniques to identify patterns. Preliminary results suggested that the data could be used as a means of exploring daily spatialtemporal behavior. The following presents insight in to the MATLAB toolboxes that can be used in detecting anomalies, that is, deviation from the normal pattern of life.

3. MATLAB toolboxes and classes for anomaly detection MATLABTM is a high-level language and interactive environment for numerical computation, visualization, and programming. MATLAB’s product family includes numerous toolboxes for parallel computing, math, statistics and optimization, control system design and analysis, signal processing and communication, image processing and computer vision, test and measurement, and other areas [7]. A key feature of these products is the ability to research a number of algorithms quickly without having to design and code the algorithm. MATLAB’s interactive environment allows end users to quickly prototype their own algo‐ rithms if existing toolboxes do not provide the desired functionality. A well-established community of users supports the exchange of algorithms for those wishing to share their research and prototypes. This research took advantage of the MATLAB community and the existing toolboxes. From the file exchange, Vladimir Bondarenko’s contribution of the class to the Twitter REST API v1.1 was used to interface with Twitter to search for and extract tweets for specific topics and geographic regions [8]. Existing MATLAB Statistics and Neural Network Toolboxes were also

5

6

MATLAB Applications for the Practical Engineer

used in this research to provide statistical algorithms and unsupervised learning methods, such as cluster analysis, for exploring data to discover hidden patterns and groupings in the data. 3.1. TWITTY — Twitter REST API interface Twitty is a useful interface that runs within MATLAB for communicating with Twitter. Methods used in Twitty are essentially wrapper functions that call the Twitter API [9]. The API caller function, callTwitterAPI(), does the main work.

Key steps toalso successfully using the Twitter APIstatistical include algorithms obtaining and Twitter credentials and using used in this research to provide unsupervised learning methods, such as clu   exploring data to(JSON) discoverparsers. hidden patterns and groupings in the data. JavaScript Object Notation The MATLAB file exchange provides a JSON parser developed by Joel Feenstra [10]. The JSON parser parses a JSON string and returns a MATLAB cell array with data. JSON objects to structures and JSON arrays are 3.1 the parsed TWITTY – TWITTER REST are APIconverted INTERFACE converted to cell arrays. Twitter credentials are easily created by registering at the Twitter site, creating an application, and retrieving consumer and access keys. These keys are required for Twitty is a useful interface that runs within MATLAB for communicating with Twitter. Methods used in Twitt running Twitty andfunctions include that specific values for:API [9]. The API caller function, callTwitterAPI(), does the main work wrapper call the Twitter

• ConsumerKey Key steps to successfully using the Twitter API include obtaining Twitter credentials and using JavaScript O

(JSON) parsers. The MATLAB file exchange provides a JSON parser developed by Joel Feenstra [10]. Th

• ConsumerSecret parses a JSON string and returns a MATLAB cell array with the parsed data. JSON objects are converted

JSON arrays are converted to cell arrays. Twitter credentials are easily created by registering at the Twitte

• AccessToken application, and retrieving consumer and access keys. These keys are required for running Twitty and incl for:

• AccessTokenSecret 

ConsumerKey

To use these keys a MATLAB script, users assign the values to their credential structure.  inConsumerSecret Then, a twitty instance can be created and methods, such as search, can be called as shown in  AccessToken  AccessTokenSecret  Figure 1. To use these keys in a MATLAB script, users assign the values to their credential structure. Then, a twitty created and methods, such as search, can be called as shown in Figure 1. % Create credentials credentials.ConsumerKey = 'YourConsumerKey’ credentials.ConsumerSecret = 'YourConsumerSecret' credentials.AccessToken = 'YourAccessToken' credentials.AccessTokenSecret = 'YourAccessTokenSecret' % Create Twitty Instance tw = twitty(credentials); % Search for World Series Related Tweets tw.search(‘World Series’);

Figure 1. Twitty search method example.

Figure 1. Twitty search method example.

Twitty provides a number of useful methods for interacting with the Twitter API. Example calls, along with a

Twitty provides a number of1.useful methods for and interacting with the Twitter API. are shown in Table Additional methods detailed descriptions are found byExample typing twitty.API at the M calls, along with a brief description, are shown in Table 1. Additional methods and detailed Method Example Call Desc descriptions are found by typing twitty.API at the MATLAB prompt. search(text string)

tw.search(‘World Series’)

search(multiple parameters)

tw.search('NFL', 'count',20, 'include_entities','true', 'geocode','39.051300,95.724660,1700mi','since_id',LastID)

Search all pu containing the Series” Search public a identifier gre LastID within Topeka, Kans the word “NFL entity data als

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

Method

Example Call

search(text string)

tw.search(‘World Series’)

Description Search all public tweets containing the words “World Series” Search public tweets having a identifier greater than LastID

search(multiple parameters)

tw.search('NFL', 'count',20, 'include_entities','true',

within 1700 miles of Topeka,

'geocode','39.051300,-95.724660,1700mi','since_id',LastI

Kansas containing the word

D)

“NFL” . Return entity data also (e.g., hashtags, URL mentions) Twit the text “Watching the

updateStatus()

tw.updateStatus(‘Watching the World Series tonight’)

World Series tonight” from your account to the public

sampleStatuses()

tw.sampleStatuses()

Return a continuous stream of random public tweets

Table 1. Twitty method examples.

3.2. MATLAB — Statistics TOOLBOXTM The Statistics Toolbox™ provides statistical and machine learning algorithms and tools for organizing, analyzing, and modeling data. Key features include data organization and management via data set arrays and categorical arrays, exploratory data analysis with the use of interactive graphics and multivariate statistics, and regression or classification for predictive modeling [11]. The toolbox also includes functions that allow users to test hypotheses more effectively by checking for autocorrelation and randomness as well as other tests, for example, t-tests, one-sample tests, and distribution tests such as Chi-square and Kolmogorov-Smirnov. Clustering algorithms available within MATLAB’s Statistical Toolbox include k-means and hierarchical approaches. Clustering algorithms are particularly useful for analyzing social networking data as they help identify natural groupings that can then be further analyzed to determine similarities or differences and make business, marketing, or other decisions. A k-means clustering algorithm forms k clusters by minimizing the mean between all cluster members. Clusters are defined by the centroid or center of each cluster. The algorithm works by moving data between clusters until the sum of the distances between each member and the centroid is minimized. MATLAB allows different distance measures to be selected. Table 2 lists the available distance measures and a brief description and example call.

7

8

MATLAB Applications for the Practical Engineer

Distance Measure

Description

Example call

Squared Euclidean

Each centroid is the mean of the

kmeans(meas,3,'dist','sqeuclidean');

points in that cluster. City Block

Sum of absolute differences.

kmeans(meas,3,'dist','cityblock');

Each centroid is the componentwise median of the points in that cluster. Cosine

One minus the cosine of the

kmeans(meas,3,'dist','cosine');

included angle between. Each centroid is the mean of the points in that cluster. Correlation

One minus the sample

kmeans(meas,3,'dist','correlation');

correlation between points. Each centroid is the component-wise mean of the points in that cluster. Hamming

Percentage of bits that differ in

kmeans(meas,3,'dist','hamming');

binary data. Each centroid is the component-wise median of points in that cluster. Table 2. MATLAB k-means clustering distance measures.

Visually displaying the results for multidimensional data can be challenging. However; the silhouette plot, available within the MATLAB Statistics Toolbox, displays a measure between 0 and 1, representing how close each point in one cluster is to the points in the neighboring clusters. Values close to 1 indicate points are distant from neighboring clusters whereas values close to 0 indicate points are not distinctly different from one cluster or another. Negative values indicate points that are most likely assigned to the wrong cluster. In the following example, cidx2 represents the cluster index for each given sample of data. Assuming cidx2 was returned from the call to the k-means function, an example call of the silhouette function is: silhouette(meas,cidx2,'sqeuclidean');

Hierarchical clustering groups data to create a tree structure consisting of multiple levels. Users can prune parts of the tree depending upon the application and level of detail required. The links between data are represented as upside-down U-shaped lines with the height of each line indicating the distance between the data. This height is known as the cophenetic correla‐ tion distance between the two objects. MATLAB has a function that measures this distance.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

9

The closer the value of the cophenetic correlation coefficient is to 1, the more accurately the clustering solution reflects your data.

Distance measures, similar to the k-means clustering, need to be selected to generate the linkages of the tree. Figure 2 shows the code to create a tree based on the Euclidean distance for a matrix named ‘meas’ and displays the results in clustering, a tree andneed thentocalculates cophenetic Distance measures, similar to the k-means be selectedthe to generate the linkages of the tre shows the code to create a tree based on the Euclidean distance for a matrix named ‘meas’ and displays the correlation. and then calculates the cophenetic correlation.

eucD = pdist(meas,'euclidean'); clustTreeEuc = linkage(eucD,'average'); [h,nodes] = dendrogram(clustTreeEuc,0); measCophenet = cophenet(clustTreeEuc,eucD);

Figure 2. MATLAB code to calculate cophenetic correlation.

Figure 2. MATLAB code to calculate cophenetic correlation.

3.3

MATLAB – NEURAL NETWORK TOOLBOXTM

TM functions and apps for modeling complex nonlinear systems that are n Neural Network Toolbox provides 3.3. MATLAB – The Neural Network TOOLBOX

modeled with a closed-form equation. Primary features include data fitting, clustering, and pattern recognition future events, both supervised and unsupervised network architectures, and training algorithms, such as gra

The Neural Network ToolboxTM provides functions and apps for modeling complex nonlinear and conjugate gradient methods, to help automatically adjust the network’s weights and biases [12]. Also inc systems that arepreprocessing not easily modeled with a closed-form includetraining and enab and post-processing functions thatequation. improve thePrimary efficiencyfeatures of neural network analysis and of network performance by reducing the dimensions of theboth input supervised vectors usingand principal component data fitting, clustering, pattern recognition to forecast future events, normalizing the mean and standard deviation of the training set. unsupervised network architectures, and training algorithms, such as gradient descent and MATLAB’s Neural Toolbox provides self-organizing maps for both unsupervised and supervised clu conjugate gradient methods, to Network help automatically adjust the network’s weights and biases organizing maps retain topological information for similar classes and provide reasonable classifiers. Self-org [12]. Also included are preprocessing and post-processing functions that improve the efficien‐ can be used for multidimensional data with complex features, making them attractive for analysis of social ne cy of neural network training and detailed analysiscan ofbenetwork performance by on the Twitter-ex The self-organizing maps enable functionality within MATLAB used to cluster tweets based reducing the dimensions of the input vectors using principal component analysis and nor‐ Self-organizing maps learn to classify input vectors according to how they are grouped in the input space. Th malizing the mean and standard of the training set. topologies. Within MATLAB, topologies for gridtop, h self-organizing map deviation can be arranged based on specific

random are possible. Results are not affected by the choice of topology used; however, visual inspection of t

MATLAB’s Neural provides grid self-organizing both has unsupervised gridtopfor topology neurons evenly spaced in ma best Network achieved byToolbox use of a hexagonal (hextop) [13]. Amaps example, an 80-neuron map may be arranged for in ansimilar 8 × 10 gridtop topology. and supervised dimensions. clustering.For Self-organizing maps self-organizing retain topological information topology is similar, but the shape of the grid can best be described as a hexagon. A random topology has the classes and provide reasonable classifiers. Self-organizing maps can be used for multidimen‐ randomly located. Figure 3 illustrates gridtop, hextop and random topologies for a 3 × 4 set of neurons. sional data with complex features, making them attractive for analysis of social networking data. The self-organizing maps functionality within MATLAB can be used to cluster tweets based on the Twitter-extracted fields. Self-organizing maps learn to classify input vectors according to how they are grouped in the input space. The neurons in a self-organizing map can be arranged based on specific topologies. Within MATLAB, topologies for gridtop, hextop and random are possible. Results are not affected by the choice of topology used; however, visual inspection of the data space is best achieved by use of a hexagonal grid (hextop) [13]. A gridtop topology has neurons evenly spaced in matrix of specific dimensions. For example, an 80-neuron self-organizing map may be arranged in an 8 × 10 gridtop topology. A hextop topology is similar, but the shape of the grid can best be described as a hexagon. A random topology has the neurons randomly located. Figure 3 illustrates gridtop, hextop and random topologies for a 3 × 4 set of neurons.

10

MATLAB Applications for the Practical Engineer

Figure 3. Gridtop, hextop, and randtop self-organizing map topologies.

The self-organization map routine can be run in MATLAB by calling the selforgmap function. Once a map is set up, it can be trained using specific input, and results can be displayed Figure 3. Gridtop, hextop, and randtop map topologies.map, visually. Figure 4 shows a code example to generate a 3 self-organizing × 4 neuron self-organizing conduct training, and display the results.

The self-organization map routine can be run in MATLAB by calling the selforgmap function. trained using specific input, and results can be displayed visually. Figure 4 shows a code ex self-organizing map, conduct training, and display the results. net = selforgmap([3 4]); net = train(net,meas); view(net); y = net(meas); classes = vec2ind(y);

 

Figure 4. MATLAB code example to generate a 3 × 4 neuron self-organizing map, conduct tr

Figure 4. MATLAB code example to generate a 3 × 4 neuron self-organizing map, conduct training, and display the results.

The Neural Network Toolbox also contains the nnstart tool, which provides an interactive too algorithms and training iterations. As shown in Figure 5, the interactive tool allows the user t The Neural Network Toolbox contains the nnstart tool, which provides an interactive tool train for also a specific number of epochs.

for loading data and selecting algorithms and training iterations. As shown in Figure 5, the interactive tool allows the user to select run options and then train for a specific number of epochs.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

Figure 5. MATLAB’s nnstart interactive interface.

The addition of the MATLAB’s Parallel Computing Toolbox™ allows the user to speed up training and handle large data sets by distributing computations and data across multiple processors and graphic processing units (GPUs). The ability to spin up multiple central processing units (CPUs) and GPUs from cloud services such as AWS makes this toolbox attractive. The user may be limited by the availability of MATLAB licenses for this option. MATLAB and AWS are currently working out a process to make this more feasible with a payas-you-go license feature that could be quite attractive for businesses that do not have the resources to purchase hardware and software. 3.4. Deploying MATLAB on AWS Cloud service providers such as AWS can significantly reduce the startup and maintenance costs for high-performance computers needed for efficient running of complex math, statistics, and optimization algorithms. AWS has free and pay tier options that will fit most budgets. Previous research has shown the cost of running multiple high-performance servers for many

11

12

MATLAB Applications for the Practical Engineer

hours was just a few dollars [14]. As long as users have a valid MATLAB license, they can install MATLAB on the AWS machine(s) and run the data and analysis. Once the runs are complete, users can disable the MATLAB license on the AWS machine. If users need to run the experiment again, they can spin up the server and enable the license to continue the experiment.

4. Experiment The experiment consisted of running MATLAB code on commercial cloud architecture (AWS) to collect publically available tweets on common keywords, such as NFL and World Series, during specific time intervals of a weekend over a large geographic region when both NFL games and World Series games were being played. Various MATLAB tools and functions were used to gather statistical information and analyze the data by use of unsupervised learning and clustering methods. For this experiment, a server running Windows2 2008 operating system was spun up in the AWS free tier. A MATLAB license was installed that included the Statistics and Neural Network toolboxes. The twitty.m and parse_json.m files were also uploaded to allow MATLAB to call the Twitter API. After successful installation, a MATLAB script file (m file) was created to provide the proper Twitter credentials, and then search for tweets related to two different sports events, the World Series and NFL games. The specific search strings used included “NFL” and “World Series”. Four different 3-hour data collections were made over a period of two days. Collection times included Saturday from 5 PM to 8 PM and 8 PM to 11 PM and Sunday from 1 PM to 4 PM and 5 PM to 8 PM during the weekend of the 2013 World Series. During this time frame, a World Series event took place on both Saturday and Sunday evenings, and a series of NFL games occurred on Sunday afternoon and evening. A large geographic collection area, centered over a 1700 mile radius about Topeka, Kansas, was used for this experiment. This geographic radius included most of the continental United States and parts of Canada and Mexico. The search parameters included the Twitter entities for storing User-, URL-, and hashtag-mentions. Sampling occurred every 60 seconds over the 3-hour window with Twitter ID information being used to eliminate duplicate tweets in the time frame. A sample call for the search is as follows: S=tw.search('NFL', 'count',20, 'include_entities','true','geocode','39.051300,-95.724660,1700mi', 'since_id',LastID);

In this call, Twitter is being searched for up to 20 tweets that include the text “NFL” that originated within 1700 miles of Topeka, Kansas. The LastID variable was updated after each 2 Registered trademark of Microsoft Corporation.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

iteration to gather only recent tweets beyond a specific ID. Appendix A shows the MATLAB m file that was run for the experiment. Statistics and analysis were run at the end of the data collection using existing MATLAB tools and functions provided in the Statistics and Neural Network Toolboxes. Within the Statistics Toolbox, the k-means and hierarchical clustering approaches and associated visual displays were used. The distance measure used for all k-means calculations was the default selection of Squared Euclidian. An example call used for three clusters with Squared Euclidian distance is as follows: kmeans(meas,3,'dist','sqeuclidean');

The self-organizing map functions within the Neural Network Toolbox were used to cluster the Twitter collected data. In addition to the existing Statistical and Neural Network functions, some simple data manipulation algorithms were used to extract the time between tweets. As shown in Figure 6, the following code example determines and stores the time between tweets found in MATLAB cell array named DatesQ2.

for i=1:length(DatesQ2)-1 time1 = [str2num(DatesQ2{i}(26:30)) 10 str2num(DatesQ2{i}(9:10)) str2num(DatesQ2{i}(12:13)) str2num(DatesQ2{i}(15:16)) str2num(DatesQ2{i}(18:19))]; time2 = [str2num(DatesQ2{i+1}(26:30)) 10 str2num(DatesQ2{i+1}(9:10)) str2num(DatesQ2{i+1}(12:13)) str2num(DatesQ2{i+1}(15:16)) str2num(DatesQ2{i+1}(18:19))]; deltaWorldSeries(i) = abs(etime(time1,time2)); end

Figure 6. MATLAB code to determine time between tweets.

Figure 6. MATLAB code to determine time between tweets.

In this code, time {i} and time {i+1} values are constructed using the default time structure of the Twitter collection and converted to a MATLAB time structure. The dates are then differenced using MATLAB’s built-in elapsed time function

In this code, time andcontinues time {i+1} values are constructed using the default time structure of (etime). This{i} process to calculate the time delta for all consecutive tweets for all search queries. the TwitterAppendix collection converted a MATLAB time structure. The dates are then differ‐ B showsand the Data analysis m fileto used for this experiment.   enced using MATLAB’s built-in elapsed time function (etime). This process continues to 5. INITIAL RESULTS AND DISCUSSION calculate the time delta for all consecutive tweets for all search queries. Results were provided using MATLAB’s visual and plotting functions. To better understand the frequency, or popularity, of the two search terms for each time period during the experiment, some baseline descriptive statistics were calculated.

Appendix B shows the Data analysis m file used for this experiment. 5.1

BOXPLOTS

Boxplots are a convenient way for visualizing the interquartile range, average and outlier data. Figures 7–10 represent boxplots for the 3-hour experiment windows from Saturday 5 PM to 8 PM, Saturday 8 PM to 11 PM , Sunday 1 PM to 4 PM, and Sunday 5 PM to 8 PM, respectively.

5. Initial results and discussion

Results were provided using MATLAB’s visual and plotting functions. To better understand the frequency, or popularity, of the two search terms for each time period during the experi‐ ment, some baseline descriptive statistics were calculated.

 

13

14

MATLAB Applications for the Practical Engineer

5.1. Boxplots Boxplots are a convenient way for visualizing the interquartile range, average and outlier data. Figures 7–10 represent boxplots for the 3-hour experiment windows from Saturday 5 PM to 8 PM, Saturday 8 PM to 11 PM, Sunday 1 PM to 4 PM, and Sunday 5 PM to 8 PM, respectively.

Figure 7. Boxplot for NFL versus World Series time between tweets, Saturday 5 PM to 8 PM.

Figure 8. Boxplot for NFL versus World Series time between tweets, Saturday 8 PM to 11 PM.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

Figure 9. Boxplot for NFL versus World Series time between tweets, Sunday 1 PM to 4 PM.

Figure 10. Boxplot for NFL versus World Series time between tweets, Sunday, 5 PM to 8 PM.

Reviewing the boxplots in Figures 7–10 reveals a couple of interesting trends. In most cases, outliers outside of the interquartile ranges exist; however, as a sporting event gets closer to start time, the frequency of tweets increase. Note the narrowing of the boxes for both the NFL

15

16

MATLAB Applications for the Practical Engineer

and World Series plots during or approaching the actual game events. By Sunday late afternoon, there was significant activity in both the NFL and World Series tweets in terms of frequency. The descriptive statistics were very similar, making it difficult to see any difference, in terms of tweet frequency between NFL and World Series fans. 5.2. XY plots MATLAB was also used to plot the time between tweets for each collection period and search query. This view, although somewhat cluttered, can be used to quickly compare the number of tweets and any other patterns obvious in an xy plot. The time between tweets is plotted in Figures 11–14, for each of the four collection periods. The most revealing information from these plots is the noticeable increase in the number of tweets as the events draw closer. On Saturday, the number of tweets related to the World Series was higher than NFL tweets. This seems reasonable since a World Series game took place on Saturday but no NFL games occurred until Sunday. On Sunday, with multiple NFL games taking place, the number of NFL tweets was larger than World Series-related tweets. However, the number and frequency of World Series tweets was also high due to the Sunday evening World Series event, but not as significant as the NFLrelated tweet activity.

Figure 11. XY plot for NFL versus World Series time between tweets, Saturday 5 PM to 8 PM.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

Figure 12. XY plot for NFL versus World Series time between tweets, Saturday 8 PM to 11 PM.

Figure 13. XY plot for NFL versus World Series time between tweets, Sunday 1 PM to 4PM.

17

18

MATLAB Applications for the Practical Engineer

Figure 14. XY plot for NFL versus World Series Time between tweets, Sunday 5 PM to 8 PM.

5.3. Histograms The histogram functionality within MATLAB was used to get a better idea of the quantity and length of the time between tweets for each of the collection periods. As shown in Figures 15 and 16, on Saturday, the number of tweets with time difference values less than 10 seconds was much greater for World Series than NFL tweets. In all histograms, a handful of tweets had a time between tweets of greater than 100 seconds. The shift to NFL interest is very visible on Sunday in the histograms shown in Figures 17 and 18. Note the significant increase in the both the number and quantity of tweets occurring in less than 10 second intervals for NFL tweets. However, in Figure 18, the histograms for Sunday late afternoon and evening reveal similar results for the time between tweets for both NFL and World Series texts.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

Figure 15. Histogram plot for NFL versus World Series time between tweets, Saturday 5 PM to 8 PM.

Figure 16. Histogram plot for NFL versus World Series time between tweets, Saturday 8 PM to 11 PM.

19

20

MATLAB Applications for the Practical Engineer

Figure 17. Histogram plot for NFL versus World Series time between tweets, Sunday 1 PM to 4 PM.

Figure 18. Histogram plot for NFL versus World Series time between tweets, Sunday 5 PM to 8 PM.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

5.4. Scatter plots Twitter queries provide additional data beyond the time and text of the tweet. For example, the number of friends, number of followers, user screen name, date of account creation, and time zone of user are also available and easily extract from the queries. This information can be used to determine if any significant differences or similarities exist among users of Twitter. Within MATLAB, scatter plots can be used to visually identify clusters of data. Figures 19–22 illustrate the use of scatter plots to compare the number of friends and number of followers for NFL (blue Xs) and World Series (red circles) tweeters during each of the four data collection periods. Although a significant amount of overlap exists in the groups, particularly for small numbers of friend and followers, outliers can be identified in all four plots. Additional variables and analysis may be needed to further isolate these groups.

Figure 19. Scatter plot for friends versus followers, Saturday 5 PM to 8 PM.

21

22

MATLAB Applications for the Practical Engineer

Figure 20. Scatter plot for friends versus followers, Saturday 8 PM to 11 PM.

Figure 21. Scatter plot for friends versus followers, Sunday 1 PM to 4 PM.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

Figure 22. Scatter plot for friends versus followers, Sunday 5 PM to 8 PM.

5.5. Silhouette plots Running the k-means algorithm in MATLAB for the Twitter feature sets further reveals the similarity in both NFL and World Series sets of tweets for all four time collection periods. Figures 23–26 illustrate the results of the applying the k-means algorithm for three clusters. In all cases, one very large cluster with silhouette values very close to 1 is observed. However, the other two clusters are very small with both negative and lower positive silhouette values. This indicates that the separation into unique clusters is difficult for this set of data and features. Additional selection of features and analysis would be needed to better identify clusters of similar Twitter users.

23

24

MATLAB Applications for the Practical Engineer

Figure 23. Silhouette plot of NFL versus World Series tweet features, Saturday 5 PM to 8PM.

Figure 24. Silhouette plot of NFL versus World Series tweet features, Saturday 8 PM to 11 PM.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

Figure 25. Silhouette plot of NFL versus World Series tweet features, Sunday 1 PM to 4 PM.

Figure 26. Silhouette plot of NFL versus World Series tweet features, Sunday 5 PM to 8 PM.

25

26

MATLAB Applications for the Practical Engineer

5.6. Hierarchical plots Hierarchical clustering provides an additional visualization of possible cluster separation for the tweet data. In all cases, the number of cluster and representatives within that cluster matches the k-means clustering results. Specifically, one large cluster was identified along with at most one or two very small clusters as illustrated in Figures 27–30 where the left portion of the tree has many nodes whereas the right portion just has one or two. The size of the cluster grows over the collection periods with the final collection period on Sunday resulting in the largest cluster with the least amount of additional clusters. For example, on Sunday, the tree’s one small branch is evident on the right side of the graph with the larger branch on the left side containing most observations.

Figure 27. Hierarchical clustering of NFL versus World Series tweet features, Saturday 5 PM to 8 PM.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

Figure 28. Hierarchical clustering of NFL versus World Series tweet features, Saturday 8 PM to 11 PM.

Figure 29. Hierarchical clustering of NFL versus World Series tweet features, Sunday 1 PM to 4 PM.

27

28

MATLAB Applications for the Practical Engineer

Figure 30. Hierarchical clustering of NFL versus World Series tweet features, Sunday 5 PM to 8 PM.

5.7. Self-Organizing Maps (SOM) plots MATLAB’s Neural Network Toolbox was also used to cluster the tweets through the selforganizing maps function. A hextop topology of 10 ×10 nodes was selected with 200 training epochs. After training, several visualizations are available to help determine how the input is distributed across the nodes. The locations of the data points and the weight vectors are shown by selecting the weight positions plot. With this display, only two weights can be shown at one time. Figure 31 shows that most of the data points cluster in one area and they are not very well distributed. This type of clustering was observed in both the hierarchical and k-means results.

Figure 31. SOM weight positions, Saturday 5 PM to 8 PM.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

5.7.1. SOM weight distance plots The SOM neighbor weight distances plot, using multiple dimensions, provides more infor‐ mation. In Figures 32–35 SOM neighbor weight distance plots are displayed for each of the four time collections. The following diagram colors and description should be used when interpreting these plots: 1.

Neurons are represented by blue hexagons

2.

Red lines connect neighboring neurons

3.

Dark-colored regions represent larger distances between neurons

4.

Light-colored regions represent smaller distances between neurons

Figures 32–35 show one large cluster, with small distances between the member records, present in all four collection periods. One or two very small clusters are also present, with relatively large distances between the member records.

Figure 32. SOM weight distances, Saturday 5 PM to 8 PM.

29

30

MATLAB Applications for the Practical Engineer

Figure 33. SOM weight distances, Saturday 8 PM to 11 PM.

Figure 34. SOM weight distances, Sunday 1 PM to 4 PM.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

Figure 35. SOM weight distances, Sunday 5 PM to 8 PM.

5.7.2. SOM weight plane plots The SOM weight plane plots are used to visualize the strength of weights that connect each input to each of the neurons. For our experiment, five inputs were used; therefore, five subplots were generated for each experiment. The five input features included, for each tweet, the number of user mentions, URL mentions, hashtag mentions, followers, and friends. Lighter colors in the plots represent larger weights whereas darker colors represent smaller weights. Similar connection patterns of the inputs indicate a high correlation. Weight plane plots are shown for each of the collection times in Figures 36–39. Inputs 4 and 5 appeared to be similar in all collection times and were interpreted as highly correlated. Input from variables 1, 2, and 3 seemed to contribute the smallest amount of cluster separation in the data sets as they appear to be the least similar and less correlated. This seems reasonable because the number of friends and followers do seem to be correlated with Twitter users as a large number of friends also have a large number of followers. The information from the user-, URL-and hashtag-mentions shows some promise as the maps show these as not being highly correlated. Additional features and analysis are recommended to enhance the differences in these maps and to perform better clustering.

31

32

MATLAB Applications for the Practical Engineer

Figure 36. SOM weight planes, Saturday 5 PM to 8 PM.

Figure 37. SOM weight planes, Saturday 8 PM to 11 PM.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

Figure 38. SOM weight planes, Sunday 1 PM to 4 PM.

Figure 39. SOM weight planes, Sunday 5 PM to 8 PM.

33

34

MATLAB Applications for the Practical Engineer

5.7.3. SOM sample hits plots One additional and useful visual plot provided by the MATLAB SOM function is the SOM sample hits plot. The sample hits plot counts the number of data records associated with each neuron. In an ideal situation, a relatively even distribution across the neurons is desired. Figures 40–43 show the distribution across the neurons was not even in our experiment. In most cases, the distribution was clustered in one area, which indicates very similar data without much separation. Even though the location of concentration seems to shift from one area of the map to another for each time frame, the results are essentially the same for each collection time. The exception is the increase in concentration of data around the “heavyhitting” neurons on Sunday as both NFL and World Series events were scheduled.

Figure 40. SOM sample hits, Saturday 5 PM to 8 PM.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

Figure 41. SOM sample hits, Saturday 8 PM to 11 PM.

Figure 42. SOM sample hits, Sunday 1 PM to 4 PM.

35

36

MATLAB Applications for the Practical Engineer

Figure 43. SOM sample hits, Sunday 5 PM to 8 PM.

6. Conclusion This research used MATLAB tools to extract and analyze social networking data sets and leverage cloud technologies and infrastructures. AWS was used to spin up a Windows 2008 server and install MATLAB and its associated toolboxes for data mining and statistics. In addition, a community MATLAB m file interfaced with the Twitter API to search specific text queries and retrieve associated data. The setup process on AWS was straightforward and provided a cost-effective and free solution for the hardware and operating system. The MATLAB license was still needed in this implementation to be able to use the toolboxes and associated m files.

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

AWS provided cost savings, including all server-related hardware and software. Since MATLAB was used, the development costs for some very specialized analysis visualization tools were also reduced. Further cost savings are envisioned as MATLAB provides cloud options for running its software for short durations using a pay-as-you-go cost model. At the time of publishing, access to MATLAB Distributed Computing Server on the cloud is available as part of an Early Adopter Program for MATLAB Distributed Computing Server on Elastic Compute Cloud (EC2). Future efforts could take advantage of this option to speed up implementations, run algorithms in parallel, and possibly further reduce licensing costs. The twitty.m and parse_json.m community files were successfully used to interface with the Twitter API to search for tweets related to specific queries including “NFL” and “World Series” sporting events. In addition to the tweet texts, other information used in the experiment included user data such as number of friends, number of followers, time between tweets, number of URL mentions, number of Hashtag mentions, and number of user mentions in each text. The MATLAB Statistics and Neural Network Toolboxes were then used to extract and display descriptive statistics and perform unsupervised clustering using k-means, hierarchical, and self-organized maps. Initial analysis of the visualization and data output revealed that sports events and the associated popularity and frequency of tweets increases as the time of the event gets closer. The frequency of tweets also persisted throughout the event. Tweet users were also shown to have similar characteristics and profiles and would be difficult to separate based on just a few features. Referencing the SOM weight distance plots, smaller clusters were observed in all the plots, however, the distance among the members of the cluster were large compared to the one large cluster observed in each plot. Additional research with larger data sets and more robust features are needed to form predicative models.

7. Recommendations This research demonstrated the use of MATLAB for social networking site analysis using AWS servers to reduce costs is feasible, cost-effective, and efficient. Further research is recommend‐ ed to Investigate MATLAB’s AWS Parallel computing license options and costs. It is believed even more cost-savings could be realized with a pay-as-you-go model. This is particularly attractive for smaller companies and start-ups that might have limited financial resources, yet have the personnel skills to conduct some excellent research. This initial experiment collected data from four different, but relatively short, time periods for two sports-related tweets. Only a handful of features were used to identify clusters and perform statistical analysis. We recommend expanding the time frame, number of tweet queries, and features to be able to provide further insight and understanding from the information available from social networking sites, such as Twitter.

37

38

MATLAB Applications for the Practical Engineer

Appendix A: Example MATLAB Code % Test Twitty call account % Clear variables clear; % Set up the credentials based on Twitter account credentials.ConsumerKey = 'YourKey'; credentials.ConsumerSecret = 'Yoursecret'; credentials.AccessToken = yourAccess; credentials.AccessTokenSecret = 'youraccesssecret'; % Create Tweet instance based on credentials % tw = twitty(credentials); % Create a Tweet example % Pick Topeka Kansas for Center of Continental U.S. % Set the last ID to Any value to start LastID = 123456; LastIDQ2 = 123456; NumIts = 180; TweetCnt=0; % Initialize to 0 counts TweetCntQ2=0; PauseTime = 60; % Start a timer Tic; % Loop several times gathering data for j=1:NumIts j % clear old data clear S; clear Q2; % Time out issues establish new tw each iteration tw = twitty(credentials); % Run the First Query S = tw.search('NFL', 'count',20, 'include_entities','true','geocode','39.051300,-95.724660,1700mi','since_id',LastID); % Let us see how many we got StatCnt = length(cellfun('ndims',S{1,1}.statuses)); for i=1:StatCnt TweetCnt = TweetCnt+1; MyTweets=S{1,1}.statuses{1,i}.text ; % Put the data in Tweets{TweetCnt}=MyTweets; % Put the dates/time in Dates{TweetCnt} = S{1,1}.statuses{1,i}.created_at; % Gather other stuff FollowerCount(TweetCnt) = S{1,1}.statuses{1,i}.user.followers_count; FriendsCount(TweetCnt) = S{1,1}.statuses{1,i}.user.friends_count; TZones = S{1,1}.statuses{1,i}.user.time_zone; TimeZones{TweetCnt} = TZones; % Get screen name SName = S{1,1}.statuses{1,i}.user.screen_name;

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

ScreenName{TweetCnt} = SName; % Get the User Mentions Count UserMentionCnt{TweetCnt} = length(cellfun('ndims',S{1,1}.statuses{1,i}.entities.user_mentions)); URLMentionCnt{TweetCnt} = length(cellfun('ndims',S{1,1}.statuses{1,i}.entities.urls)); HashTagsMentionCnt{TweetCnt} = length(cellfun('ndims',S{1,1}.statuses{1,i}.entities.hashtags)); end % Get the lastID to eliminate Dups if (StatCnt > 0) LastID = S{1,1}.statuses{1,1}.id_str; end % Run the World Series Query Q2 = tw.search('World Series', 'count',20, 'include_entities','true','geocode','39.051300,-95.724660,1700mi','since_id',LastIDQ2); % Let us see how many we got StatCntQ2 = length(cellfun('ndims',Q2{1,1}.statuses)); for i=1:StatCntQ2 TweetCntQ2 = TweetCntQ2+1; MyTweetsQ2=Q2{1,1}.statuses{1,i}.text ; % Put the data in TweetsQ2{TweetCntQ2}=MyTweetsQ2; % Put the dates/time in DatesQ2{TweetCntQ2} = Q2{1,1}.statuses{1,i}.created_at; % Gather other stuff FollowerCountQ2(TweetCntQ2) = Q2{1,1}.statuses{1,i}.user.followers_count; FriendsCountQ2(TweetCntQ2) = Q2{1,1}.statuses{1,i}.user.friends_count; TZonesQ2 = Q2{1,1}.statuses{1,i}.user.time_zone; TimeZonesQ2{TweetCntQ2} = TZonesQ2; % Get screen name SNameQ2 = Q2{1,1}.statuses{1,i}.user.screen_name; ScreenNameQ2{TweetCntQ2} = SNameQ2; % Get the User Mentions Count UserMentionCntQ2{TweetCntQ2} = length(cellfun('ndims',Q2{1,1}.statuses{1,i}.entities.user_mentions)); URLMentionCntQ2{TweetCntQ2} = length(cellfun('ndims',Q2{1,1}.statuses{1,i}.entities.urls)); HashTagsMentionCntQ2{TweetCntQ2} = length(cellfun('ndims',Q2{1,1}.statuses{1,i}.entities.hashtags)); end % Get the lastID to eliminate Dups if (StatCntQ2 > 0) LastIDQ2 = Q2{1,1}.statuses{1,1}.id_str; end % Pause a few seconds pause(PauseTime); end

39

40

MATLAB Applications for the Practical Engineer

Appendix B: Example Data Analysis MATLAB M-File % Combine the Data into a vector measWS = [cell2mat(UserMentionCntQ2)' cell2mat(URLMentionCntQ2)' cell2mat(HashTagsMentionCntQ2)' FollowerCountQ2' FriendsCountQ2']; measNFL= [cell2mat(UserMentionCnt)' cell2mat(URLMentionCnt)' cell2mat(HashTagsMentionCnt)' FollowerCount' FriendsCount']; % Combine meas = [measWS;measNFL]; % Now let's do some Stats and Clustering [cidx2,cmeans2] = kmeans(meas,2,'dist','sqeuclidean'); figure; [silh2,h] = silhouette(meas,cidx2,'sqeuclidean'); % Look for maximum number of clusters [cidx3,cmeans3] = kmeans(meas,3,'display','iter'); % Now do some clustering an figure; [silh3,h] = silhouette(meas,cidx3,'sqeuclidean'); % Some plots Might want to plot with NFL versus World Series % Hierarchical Clustering eucD = pdist(meas,'euclidean'); clustTreeEuc = linkage(eucD,'average'); myCop = cophenet(clustTreeEuc,eucD) figure; [h,nodes] = dendrogram(clustTreeEuc,0); set(gca,'TickDir','out','TickLength',[.002 0],'XTickLabel',[]); % Reduce nodes figure; [h,nodes] = dendrogram(clustTreeEuc,12); toc % World Series tweet Gaps for i=1:length(Dates)-1 time1 = [str2num(Dates{i}(26:30)) 10 str2num(Dates{i}(9:10)) str2num(Dates{i} (12:13)) str2num(Dates{i}(15:16)) str2num(Dates{i}(18:19))]; time2 = [str2num(Dates{i+1}(26:30)) 10 str2num(Dates{i+1}(9:10)) str2num(Dates{i+1}(12:13)) str2num(Dates{i+1}(15:16)) str2num(Dates{i+1} (18:19))]; deltaWorldSeries(i) = abs(etime(time1,time2)); end % NFLtweet gaps for i=1:length(DatesQ2)-1 time1 = [str2num(DatesQ2{i}(26:30)) 10 str2num(DatesQ2{i}(9:10)) str2num(DatesQ2{i}(12:13)) str2num(DatesQ2{i}(15:16)) str2num(DatesQ2{i}(18:19))]; time2 = [str2num(DatesQ2{i+1}(26:30)) 10 str2num(DatesQ2{i+1}(9:10)) str2num(DatesQ2{i+1}(12:13)) str2num(DatesQ2{i+1}(15:16)) str2num(DatesQ2{i+1} (18:19))]; deltaNFL(i) = abs(etime(time1,time2)); end figure (1) subplot(2,1,1) plot (deltaNFL,'r');

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

ylabel('Time between Tweets') xlabel ('Tweet Number') title('"NFL" Time between Tweets'); grid; subplot(2,1,2) plot(deltaWorldSeries, 'b'); ylabel('Time between Tweets') xlabel ('Tweet Number') title('"World Series" Time between Tweets'); grid; % Stats % Quick histograms for the 2 delta times xvalues = -1:2:800; figure (2) subplot (1,2,1) hist(deltaNFL,xvalues); grid; title ('Time between "NFL" Tweets') axis ([-1 200 0 200]); subplot (1,2,2) hist(deltaWorldSeries,xvalues); grid; title ('Time between "World Series" Tweets') axis ([-1 200 0 200]); % Box plots for Descriptive visualization figure(3) subplot(1,2,1) boxplot(deltaNFL) title ('Time between "NFL" Tweets') grid; subplot(1,2,2) boxplot(deltaWorldSeries) title ('Time between "World Series" Tweets') grid; % Group Scatter % Put everything together MyGroup1 = zeros(1,length(HashTagsMentionCnt)); MyGroup2 = ones(1,length(HashTagsMentionCntQ2)); MyGroups = [MyGroup1 MyGroup2]; A1 = cell2mat(HashTagsMentionCnt); B1 = FollowerCount; C1 = FriendsCount; D1 = cell2mat(UserMentionCnt); E1 = cell2mat(URLMentionCnt); A2 = cell2mat(HashTagsMentionCntQ2); B2 = FollowerCountQ2; C2 = FriendsCountQ2; D2 = cell2mat(UserMentionCntQ2); E2 = cell2mat(URLMentionCntQ2); MyVars = [A1' B1' C1' D1' E1'; A2' B2' C2' D2' E2']; varNames = {'HashTagMention' 'FollowerCount' 'FriendCount' 'URLMention' 'HashTagMention'};

41

42

MATLAB Applications for the Practical Engineer

figure(4) gplotmatrix(MyVars(:,2), MyVars(:,3),MyGroups,['r' [],'off','hist',['Followers'],['Friends']); grid;

'b'],['O'

'X'],

Author details Kelly Bennett1 and James Robertson2 1 U.S. Army Research Laboratory, Sensors and Electron Devices Directorate, Adelphi, MD, USA 2 Clearhaven Technologies LLC, Severna Park, MD, USA

References [1] Rahman, M. Mining social data to extract intellectual knowledge. doi:10.5815/ijisa. 2012.10.02 [2] Das TK., Kumar P. BIG Data Analytics: A Framework for Unstructured Data Analy‐ sis. International Journal of Engineering Science and Technology 2013;5(2) 153-156. [3] Daehoon K., Daeyong K., Seungmin R., Eenjun H. Detecting Trend and Bursty Key‐ words Using Characteristics of Twitter Stream Data. International Journal of Smart Home 2013;7(1) 209-220. [4] Karandikar A. Clustering short status messages: A topic model based approach. Master’s thesis. University of Maryland; 2010. http://ebiquity.umbc.edu/_file_directo‐ ry_/papers/518.pdf (accessed 29 October 2013). [5] Perera R. Twitter Analytics: Architecture, Tools and Analysis, Proc. The 2010 Military Communications Conference, MILCOM2010, Oct. 31, 2010-Nov. 3, 2010, San Jose Convention Center, CA; 2010. [6] Turner and Malleson (2011). Applying geographical clustering methods to analyze geo-located open micro-blog posts: A case study of Tweets around Leeds since June 2011, https://docs.google.com/document/d/1yaRBUwJy8Cb3JWlNfJOZZRN1JnXvpdWzhRDstqwsoI/edit?pli=1#heading=h.rx5m14o2e6yw (accessed 30 September 2013). [7] MathWorks. MATLAB Products and Services. Overview. http://www.math‐ works.com/products/matlab/ (accessed 8 October 2013). [8] MathWorks. MATLAB Central. File Exchange. twitty by Vladimir Bondarenko. 30 January 2012 (Updated 12 July 2013). Interface-class to access the Twitter REST API

Knowledge Discovery and Information Extraction on the Open Internet Using MATLAB and Amazon… http://dx.doi.org/10.5772/58895

v1.1. http://www.mathworks.com/matlabcentral/fileexchange/34837-twitty/content/ twitty.m. (accessed 21 September 2013). [9] Twitter API. https://dev.twitter.com/docs/api/1.1 (accessed 12 July 2013). [10] MathWorks. MATLAB Central. File Exchange. JSON Parser. JSON Parser by Joel Feenstra. 3 July 2008 (Updated 18 June 2009). http://www.mathworks.com/matlab‐ central/fileexchange/20565-json-parser. (accessed 23 September 2013). [11] MathWorks. Statistics Toolbox. Overview. http://www.mathworks.com/products/ statistics/ (accessed 1 October 2013). [12] MathWorks. Neural Network Toolbox Overview. http://www.mathworks.com/prod‐ ucts/neural-network/ (accessed 18 September 2013). [13] Kohonen T. Self-Organizing Maps. 3rd ed. Springer Series in Information Sciences. Springer-Verlag Heidelberg: New York; 2001. [14] Bennett K, Robertson J. Signal and image processing algorithm performance in a vir‐ tual and elastic computing environment, Proc. SPIE Vol. 8734, Active and Passive Signatures IV, 87340B (2013).

43

Chapter 2

A MATLAB-based Microscope Milton P. Macedo Additional information is available at the end of the chapter http://dx.doi.org/10.5772/58532

1. Introduction When someone intends to build a laboratorial prototype of a microscope there are two major tasks. One is in the optical side which is the selection of adequate optical and mechanical components of the optical setup taking into account budget restrictions. However the major challenge is to find the adequate overall environment that enables a easy and effective integration of the different parts of the microscope in order to arrange an efficient instrument. This is a consequence of the great evolution in microscopy field in the last decades that naturally follows the progress in science and particularly in electronics and computer science. Currently the microscopes have little similarity with the general concept of an ancient microscope. In fact those former stand-alone microscopes that were used in biology laborato‐ ries at school have moved on to a complete instrumentation system. Different areas such as optics, mechanics, electronics and software have now to be integrated in order to get a digital image of an object. Ultimately a modern microscope is a user computer-controlled instrument. Surely this configures an instrumentation field where it is very attractive to use MATLAB. In this chapter it is presented a practical example of MATLAB application as the fundamental tool in a three-dimensional (3D) microscopy platform. The first stage of this research project consisted on the selection of a one-dimensional (1D) array CCD/CMOS sensor and the subsequent development of the sensor readout module. Afterwards the laboratory platform has been built. Besides the sensor readout module the main components of this benchmicroscope are the optical layout and computer software. The choice of an efficient computer software is fundamental as the configuration of acquisition parameters, control of data acquisition, visualization of microscope images and data process‐ ing are to be performed by the user in a computer. The versatility of this platform is probably the most important feature in order to accommodate such different tasks as:

48

MATLAB Applications for the Practical Engineer

• Acquisition of data from stand-alone sensor readout module-implementation of a data communication channel from sensor readout module to the computer; • Positioning control of the object stage – implementation of the interface of stage actuators to computer; • Visualization of data in real-time in a computer display; • Development of 3D reconstruction algorithms; • Development of Graphical User Interfaces (GUI) for the different applications. A brief reference has to be made to the use of the Borland C compiler as computer software development tool at an initial stage of this work. Namely when testing a different acquisition hardware using a different sensor and to test this sensor readout module. But next the integration of those different tasks enumerated above had to be accomplished. At this point MATLAB had been naturally considered owing to its versatility given by the functionalities from its core and toolboxes. In this chapter the focus is the description of MATLAB applications. However for a deeper comprehension of MATLAB functionalities implemented in these applications some details of the bench-microscope prototype have to be stated. Then firstly some hardware features as the sensor readout module and object stage positioning are reported. The particularities in image build and visualization owing to the use of a linear image sensor in this bench-microscope are also covered. On the other hand the best mode of evaluating the effectiveness of the MATLAB applications is showing the results obtained with this bench-microscope. Four applications have been developed for the implementation of image acquisition and visualization as well as for the assessment of image quality and image processing in some practical applications of this platform in materials science field. Lastly a summary of the overall functionalities of these different MATLAB applications and a discussion of the advantages of a platform that use such a diversity of integrated tools is presented.

2. Microscope implementation The challenge of assembling a microscope with the diversity of areas of knowledge that it demands had led to the development of open-source microscopy software. In this manner different research groups in universities as well as in industry work together in order to build software platforms that make easier the implementation of the different tasks in order to accomplish the acquisition of an image on a microscope. Obviously the use of an open source microscopy software should be considered whenever a new microscope setup is assembled. Amongst these open source microscopy software, MicroManager and OME (Open Microscopy Environment) are probably the better established.

A MATLAB-based Microscope http://dx.doi.org/10.5772/58532

Micro-Manager is a complete microscopy software that include control of microscope itself and associated hardware. A long list of microscope models as well as cameras, stages and other peripherals can be controlled. It runs as a plugin to ImageJ and provides an easy to use software to control microscope experiments. OME aims at the storage and manipulation of biological data. It comprises a client-server software (OMERO) for visualization, management and analysis of images and a Java library (Bio-Formats) for reading and writing biological image files. This library can be used as an ImageJ plugin, MATLAB toolbox or in our own software. ImageJ is a public domain Java image processing program. It is a very complete image tool that can be used with many image formats as well as raw-images. As this is a research project a fast and easy access to hardware in order to test other acquisition and control configurations would be important. Surely if a software platform is developed from the zero it is more versatile and flexible. On the other hand the graphical user interface (GUI) was adapted from the one used in the preliminary tests which had been developed in C/C++language. These two issues together with other particularities listed in table 1 were decisive to the choice of developing in MATLAB an entirely original and dedicated software for this project.

Micro-Manager

OME Disadvantages

Sensor and stage actuators devices not supported

Original MATLAB applications Advantages

Do not perform microscope control

Versatility and flexibility in hardware configuration

Time consuming and energy

Optimized to work with biological

GUI layout previously developed in

investment in Micro-Manager

data

C/C++ language

No previous experience in using

No previous experience in using

Experience in programming

ImageJ

ImageJ

languages / MATLAB

Table 1. Advantages and disadvantages of using an open source microscopy software or develop original Matlab applications for this project.

3. Description of bench-microscope The purpose of this work was to build a platform to develop and test algorithms able to obtain 3D images. The own microscope optics is based on a linear-array image sensor. At an initial stage an hardware previously developed in our research team, named PAF (Photodiode Array Fluorometer), had been used. After the implementation of a few and plain adjustments in the PAF software the first tests were performed. A scheme of this system is shown in figure 1.

49

applications for this project.

3. Description of bench-microscope 50

MATLAB Applications for the Practical The purpose of this work was to buildEngineer a platform to develop and test algorithms able to obtain 3D images. The own microscope optics is based on a linear-array image sensor. At an initial stage an hardware previously developed in our research team, named PAF (Photodiode Array Fluorometer), had been used. After the implementation of a few and plain adjustments in the PAF software the first tests were performed. A scheme of this system is shown in figure 1. CCD PC PAFPC

OPTICAL SYSTEM

PAFDET

OTS

PDU

Figure 1. Test system using PAF. PAFPC e PAFDET are PAF hardware modules in computer bus and with detector, Figure 1. TestOTS system usingTranslation PAF. PAFPCStage; e PAFDET PAF hardware modules respectively; – Object PDU are – Piezoelectric Drive Unit. in computer bus and with detector, re‐ spectively; OTS – Object Translation Stage; PDU – Piezoelectric Drive Unit. Owing to the excessive pixel width in the linear-array CCD used in PAF other sensor must be searched. A new sensor in the market (LIS-1024 from Photon-Vision Systems) with 1024 pixels of 7,8 μm width had been selected. Afterwards the development of a sensor readout also enabled to assemble theCCD sensor in theinoptical blockmust diagram Owing to the excessive pixelhad width in the linear-array used PAF bench. other The sensor be showing microscope architecture is presented in figure 2 as well as a photo of the optical layout showing the sensor-

searched. A new sensor in the market (LIS-1024 from Photon-Vision Systems) with 1024 pixels of 7,8 μm width had been selected. Afterwards the development of a sensor readout had also enabled to assemble the sensor in the optical bench. The block diagram showing microscope architecture is presented in figure 2 as well as a photo of the optical layout showing the sensorreadout module and stage actuators both controlled from MATLAB. The optical layout is beyond the scope of this chapter but for clarity a brief description of sensor-readout module and positioning of the object platform is essential.

Sensor-readout module

RS 232

Sensor

Óptical Layout

PC

Three-axis Translation Stage Positioning RS 232 readout module and stage actuators both controlled from MATLAB. The optical layout is beyond the scope of this chapter but for2. clarity brief description of sensor-readout module and positioning of the object platform is essential. Figure Blocka diagram and a photo of the bench-microscope. Figure 2. Block diagram and a photo of the bench-microscope.

3.1. Sensor-readout module

3.1. Sensor-readout module

This stand-alone module is based is on based a microcontroller of PIC family (PIC16F876) fromfamily Microchip. It had been selected This stand-alone module on a microcontroller of PIC (PIC16F876) from amongst a set of similar devices as it completely meets the predefined specifications, namely: a 10-bit ADC, three timers Microchip. It had owing beentoselected a set of similar devices as it completely meets the and an high versatility an interruptamongst structure with thirteen interrupt sources. predefined namely: a 10-bit andcommunication. an high versatility owing Its weakness specifications, lies in communication options. It onlyADC, has a three USARTtimers for RS232 So sensor data to is transferred to the computer through its serial port (RS232).sources. However the optimization of the system regarding acquisition an interrupt structure with thirteen interrupt speed is not a goal of this project. Otherwise other PIC microcontroller, PIC16C745, with an USB serial port would be the

right choice. But the overall specifications of PIC16F876 are more adequate to the systems needs namely because of its Its weakness lies in communication options. It only has a USART for RS232 communication. 10-bit ADC in comparison to the 8-bit ADC in PIC16C745. So sensor data is transferred to the computer through its serial port (RS232). However the Figure 3 shows the block diagram of the module. Besides the sensor and microcontroller it contains a RS-232 driver (MAX242 from Maxim) receive/transmit signals from/to PICspeed serial port. Thisadriver put data in electric format of optimization of thethat system regarding acquisition is not goalalso of this project. Otherwise RS-232 standard and manage control signals for data communication with the computer. As serial ports have been other PIC microcontroller, PIC16C745, with an USB serial port would be the right choice. But gradually disappearing of computers in recent years there is the possibility of using an USB to RS232 adapter and Sensor readout module Sensor

RS-232 Driver

Sensor data Commands

PC

Positioning RS 232 readout module and stage actuators both controlled from MATLAB. The optical layout is beyond the scope of this chapter A MATLAB-based but for clarity a brief description of sensor-readout module and positioning of the object platform is essential. Microscope http://dx.doi.org/10.5772/58532 Figure 2. Block diagram and a photo of the bench-microscope.

the overall specifications of PIC16F876 are more adequate to the system needs namely because 3.1. Sensor-readout module ofThis itsstand-alone 10-bit ADC in comparison to the 8-bit ADC in PIC16C745. module is based on a microcontroller of PIC family (PIC16F876) from Microchip. It had been selected amongst a set of similar devices as it completely meets the predefined specifications, namely: a 10-bit ADC, three timers

Figure 3 shows the block of thewith module. Besidessources. the sensor and microcontroller it and an high versatility owing to andiagram interrupt structure thirteen interrupt contains a RS-232 (MAX242 fromIt Maxim) receive/transmit signals from/to PICdata serial Its weakness lies in driver communication options. only has that a USART for RS232 communication. So sensor is transferred the computer through port (RS232). the optimization of themanage system regarding port. This to driver also put dataitsinserial electric formatHowever of RS-232 standard and controlacquisition signals speed is not a goal of this project. Otherwise other PIC microcontroller, PIC16C745, with an USB serial port would be the for data communication with the of computer. serialadequate ports have been gradually disappearing right choice. But the overall specifications PIC16F876 As are more to the systems needs namely because of its ADC in comparison to the 8-bit ADC in PIC16C745. of10-bit computers in recent years there is the possibility of using an USB to RS232 adapter and Figure 3 this showsmodule the blockto diagram of the module.and Besides the sensor and microcontroller it contains a RS-232 driver connect a more modern powerful computer. (MAX242 from Maxim) that receive/transmit signals from/to PIC serial port. This driver also put data in electric format of RS-232 standard and manage control signals for data communication with the computer. As serial ports have been gradually disappearing of computers in recent years there is the possibility of using an USB to RS232 adapter and Sensor readout module Sensor

RS-232 Driver

Sensor data Commands

PC

Control  signals 

ADC 10-bit

3 Timers

Memory

USART

ICSP

PIC16F876

In-Circuit Serial Programming (ICSP)

connect this module to a more modern and powerful computer.

Figure sensor readout module. Figure3.3.Block Blockdiagram diagram of of the the sensor readout module. Another important functionality of this PIC microcontroller, shown in this block diagram, is called ICSP (In-Circuit Serial

Another important functionality of for thisprogramming PIC microcontroller, shown in this block diagram, is Programming). This enables an easy mode the PIC, changing its configurations, namely the sensor readout mode, without the need of take the hardware module out of optical layout. called ICSP (In-Circuit Serial Programming). This enables an easy mode for programming the PIC, changing its configurations, namely the sensor readout mode, without the need of take the hardware module out of optical layout. 3.1.1. Sensor readout The ICSP functionality also enables a fast and easy configuration of the sensor operation mode taking advantage of its flexibility that comes out from CMOS technology. It is very useful as the more adequate sensor operation mode depends on the type of application in which the bench-microscope is used. The sensor readout mode normally used was the destructive Dynamic Pixel Reset (DPR). The reset of each pixel is executed as the respective readout is concluded. This assures an equal integration time for every pixel. To ensure the correct timing of the readout start of the next set of 1024 pixels, the PIC waits by a specific control signal from the sensor. The verification of timing specifications for sensor readout is ensured by the three timers of the PIC. Therefore the internal clock of the sensor with a duty-cycle of 50%, the readout of each

51

52

MATLAB Applications for the Practical Engineer

pixel data in the second half of the respective readout time window and the acquisition time in accordance with precision specifications of the ADC are easily implemented. The acquisition cycle is based on the appropriate interrupt structure of the PIC. Figure 4 presents the timing diagram of the acquisition cycle. It is also shown the timing of data transfer to computer through RS-232 serial port. No external memory exists in this stand-alone module. As there is no way to store the data in the module memory, each pixel value of 10 bits, the result of the ADC conversion of the analog value read from each sensor pixel, has to be sent to computer till the end of the timeslot. In this case the timeslot is the time lapse from one ADC conversion to the next one. These 10 bits value from each pixel is packed in a frame with three words of 8 bits (bytes) as it is shown in figure 5. Thus each timeslot must be long enough for the USART complete the transfer of this frame. Using the Instrument Control Toolbox the configuration of the computer serial port was performed. Preliminary tests had shown that using the maximum baud rate of 59200 bps the communication errors were very scarce. In spite of this it had been considered that a baud rate of 19200 bps was the best compromise between speed and reliability. The option had been to completely avoid these error relaxing the speed goals. This lower baud rate imposes a rate of pixel sensor readout slightly above 1 kHz. This is achieved from a timeslot width (Tcycle) of nearly 900 μs. Therefore the acquisition of all 1024 sensor pixels takes around one second. In many applications of this bench microscope it is unnecessary to perform the acquisition of the 1024 pixels. Owing to the easiness of programming the PIC that arises from the ICSP functionality described above it is plain to change the sensor readout configuration namely the timing issues in order to adequate it to the specific needs of each situation. In these cases instead of wasting time for the acquisition of data without relevant information it is clearly useful the definition of a region of interest (ROI). The dimension of this ROI in image plane depends on the magnification of the objective used. One of the objective lenses in this bench microscope is 40x / 0.65 NA (numerical aperture is a number related to the width of the cone of light gathered by the lens). So the image with 1024 pixels corresponds to 200 μm in the sample (in object plane). If in one specific application this is excessive a ROI should be de‐ fined and consequently the acquisition rate is increased. For a ROI with 256 or 128 pixels the acquisition rate is increased by a factor of 4 or 8, respectively. 3.2. Positioning actuators This bench microscope is intended to be used in reflection mode. Its optical layout is an epiillumination configuration, typical of confocal microscopes. In this case light travelling from the light source to the sample has a fraction of its path in common to light reflected by the sample. Due to budget constraints it has been made the option by a stage-scanning instead of the beam-scanning configuration.

A MATLAB-based Microscope http://dx.doi.org/10.5772/58532

Tcycle

Tcycle

Timer1 CLK interrupt (sensor input) 

Tcycle Tdelay_timer1

Tdelay_timer1

Timer1 Timer0 interrupt interrupt

count 3 toTcycle 1024

count 2

count 1

CLK (sensor input) 

count 3 to 1024

count 2

count 1

Tcycle

Tcycle

Ttimer1

Tacq

Ttimer1

Tacq

Ttimer1

Tacq

Tacq

Ttimer1

Tacq

Tacq

Timer0 A/D_Go/Done interrupt bit Status

Tconv

Tconv

A/D_Go/Done A/D bitinterrupt Status

pixel 1

Vo (sensor output/ A/DRS-232 input) data link

Pixel 3 to 1024 Tconv readout

Pixel 1 and 2 Tconv readout

Tconv A/D Vo interrupt (sensor output/ A/D input)

Tconv

Pixel 3 to 1024 pixel 3 readout

pixel 2

Pixel 1 and 2 readout

pixel 1

pixel 3

pixel 2

RS-232 data link

frame 1

frame 2

the computer display. Also these positioning products use a direct drive system for a simplified mechanical design with frame 1 frame 2 no coupling, gear, belt or other expensive components. Likewise the specifications of this linear actuators in terms of resolution and repeatability comply with the purposes of the the computer display. Also these positioning products use a direct drive system for a simplified mechanical design with bench-microscope thence the by the linear actuator T-LAdata 28, with 28 mm range. Figure 4. Timing diagram of option the acquisition cycle of sensor (timea intervals not in scale). Pixel acquisition time no coupling, gear, belt or other expensive components. (T 888 μs, A/D acquisition time (Tacq) of 30 μs and conversion time (Tconv) not of 20 μs. cycle) of Figure 4. Timing diagram of the acquisition cycle of sensor data (time intervals in Pixel timeof(T cycle) Likewise the specifications of this linear actuators in terms of resolution and repeatabilityscale). comply withacquisition the purposes the ofthe 30 linear s andactuator conversion (Tconva) 28 of 20 of 888 s, A/D acquisition time (Tacq)by bench-microscope thence the option T-LAtime 28, with mms. range. Figure 4. Timing diagram of the acquisition cycle of sensor data (time intervals not in scale). Pixel acquisition time (Tcycle) Bytetime 1 (T ) of 30 s and conversion Byte 3 Bytetime 2 (Tconv) of 20 s. of 888 s, A/D acquisition acq S

O F Byte 1

9

8

7

6

5

4

3

2

9 A/D 8 7 6 5 A/D 4 Low 3 Byte 2 High Byte (8-bit) (2-bit) A/D Low Byte A/D Start Of Frame Figure 5. Format of(6-bit) the frame used in Byte RS-232 communication. High (8-bit) Figure 5. Format of the frame used(2-bit) in RS-232 communication.

StartSOf Frame O F (6-bit)

1

0

E O Byte 3

Byte 2 1

0

Not Used Not Used

E

F

End Of Frame O (6-bit) F

End Of Frame (6-bit)

The resolution (or addressability) is the distance equivalent to the smallest incremental move the device can be instructed to make. In frame other used words, it is thecommunication. linear or rotational displacement corresponding to a single microstep of Figure 5. Format of the in RS-232 movement. The resolution for T-LA products is 0.09921875 mm (or approximately 0.1 mm). motors, acousto-optic There are a wide range of positioning e.g.,incremental stepper The resolution (or addressability) is the distance devices equivalent that to theuse, smallest move the device can be The repeatability is In theother maximum the or position of the device when attempting totoreturn amicrostep position after instructed to make. words,deviation it is the in linear rotational displacement corresponding a single of deflectors (AOD), galvanometric mirrors or piezoelectric drivers. The selection oftothe actuators moving to aThe different position. Theproducts typical repeatability T-LA actuators is about 0.3 mm. Also the typical backlash, movement. resolution for T-LA is 0.09921875formm (or approximately 0.1 mm). which is the the deviation of the finalofposition that results from reversing the direction of approach, is 2.2 mm. T-LAissues: devices to control positioning the object translation stage hadwhen been based on the following The is the maximum deviation the position of the attempting to return to absolute a positionposition, after haverepeatability built in anti-backlash routines that doin not affect motion indevice the positive direction (increasing moving to a different position. The typical repeatability for T-LA actuators is about 0.3 mm. Also the typical backlash, plunger extending). For negative motion, however, the device will overshoot the desired position by roughly 600 which is theand deviation ofapproaching the final position thatthree-axis results from reversing the direction of approach, is 2.2 T-LA devices • Easiness toreturn, accommodate in requested the translation stage (Melles Griot 17mm. AMB 003); microsteps the position from below. have built in anti-backlash routines that do not affect motion in the positive direction (increasing absolute position, plunger extending). For negative motion, however, the device will overshoot the desired position by roughly 600 microsteps and return,control approaching the requested position from below. • Computer-controlled; 3.2.1. Positioning

3.2.1. Positioning • Cost effective control (nice compromise between cost and performance).

53

54

MATLAB Applications for the Practical Engineer

T-series positioning products from ZaberTM and in particular linear actuators are ready to mount in the translation stage. These computer controlled positioning products use stepper motors to achieve open loop position control. These devices turn by a constant angle called a step for every electrical impulse sent to them. This allows a system to be built without feedback, reducing total system cost. However, being incremental (as opposed to absolute) in nature, the stepper motor must initially be zeroed by going to a home sensor. As there is no encoder, the actual position of the device will become different from the position shown in the computer display. Also these positioning products use a direct drive system for a simplified mechanical design with no coupling, gear, belt or other expensive components. Likewise the specifications of this linear actuators in terms of resolution and repeatability comply with the purposes of the bench-microscope thence the option by the linear actuator TLA 28, with a 28 mm range. The resolution (or addressability) is the distance equivalent to the smallest incremental move the device can be instructed to make. In other words, it is the linear or rotational displacement corresponding to a single microstep of movement. The resolution for T-LA products is 0.09921875 mm (or approximately 0.1 mm). The repeatability is the maximum deviation in the position of the device when attempting to return to a position after moving to a different position. The typical repeatability for T-LA actuators is about 0.3 mm. Also the typical backlash, which is the deviation of the final position that results from reversing the direction of approach, is 2.2 mm. T-LA devices have built in anti-backlash routines that do not affect motion in the positive direction (increasing absolute position, plunger extending). For negative motion, however, the device will overshoot the desired position by roughly 600 microsteps and return, approaching the requested position from below. 3.2.1. Positioning control Image acquisition with this bench-microscope requires to control the positioning of the object stage. Automatic scanning in the three-axis (XYZ) is performed. One T-LA 28 actuator controlled through the RS232 port of a computer is used to implement the scanning in each axis. However the three units are connected in a daisy-chained mode thus sharing the same serial port in the computer. The configuration of the computer serial port and the control of RS232 communication was implemented in a MATLAB application. Communications settings must be: 9600 baud, no hand shaking, no parity, one stop bit. After power-up, the units in the chain will each initialize itself as unit #1 and thus each will execute the same instructions. To assign each unit a unique identifier, a renumber instruction must be issued after all the units in the chain are powered up and every time a unit is added or removed from the chain. Instructions must not be transmitted while the chain is renumbering or the renumbering routine may be corrupted. Renumbering takes less than a second, after which instructions may start to be issued over the RS232 connection.

A MATLAB-based Microscope http://dx.doi.org/10.5772/58532

All instructions consist of a group of 6 bytes. They must be transmitted with less than 10 ms between each byte. If the unit has received less than 6 bytes and then a period of more than 10 ms passes, it ignores the bytes already received. The table 2 below shows the instruction format: Byte 1

Byte 2

Byte 3

Byte 4

Byte 5

Byte 6

Unit #;

Command #

Data (least

Data

Data

Data (most

significant byte)

significant byte)

Table 2. Instruction format.

The first byte is the unit number in the chain. Unit number 1 is the unit closest to the computer, unit number 2 is next and so forth. If the number 0 is used, all the units in the chain will execute the accompanying command simultaneously. The second byte is the command number. Bytes 3, 4, 5 and 6 are data in long integer, 2’s complement format with the least significant byte transmitted first. How the data bytes are interpreted depends on the command. Most instructions cause the unit to reply with a return code. It is also a group of 6 bytes. The first byte is the device #. Byte 2 is the instruction just completed or 255 if an error occurs. Bytes 3, 4, 5 and 6 are data bytes in the same format as the instruction data byte. For some instructions in this reply it is sent the actual effective position. Therefore the communication between the computer and T-LA 28 is executed in both direc‐ tions. From the graphical user interface (GUI) of the MATLAB application the user gives an order to perform a command, e.g., reset, home, renumber, move absolute, move relative, using the instruction format presented in table 1. The slow data communication between these linear actuators and the computer, 9600 baud, is one of the most important weaknesses of these actuators. Image acquisition rate is conse‐ quently low but it is assumed that in this work the concern is not to optimize this rate. The aim of this project is to obtain microscope images from which a three-dimensional reconstruction is made. Therefore with this goal in mind the easiness in the integration of this actuators in the optical layout combined with programming versatility largely compensate the imposed low acquisition rate.

4. Image visualization Raw-data provided in one sensor readout consists of 1024 analogue values that the A/D in microcontroller converts in 10-bit digital values. And these values compose a linear image of a region-of-interest (ROI) in the object with a length that depends on the objective magnifica‐ tion. It corresponds to 200 μm and 400 μm for 40X and 20X objectives, respectively, Usually the object plane is considered as the XY plane. Hence the optical axis which is perpendicular to this plane is the z-axis. When using a linear sensor normally the x-axis

55

56

MATLAB Applications for the Practical Engineer

represents the direction of the sensor pixels. Thus the acquisition of a two-dimensional (2D) image implies the scanning of the object on the translation-stage in the other lateral direction perpendicular to sensor (y-axis). 4.1. 2D image view Two-dimensional image build is schematically shown in figure 6. This diagram illustrates its dependence on spatial sampling rate in both lateral axes. Sampling rate in x-axis is imposed by pixel width which is 7,8 μm. It is possible to join contiguous pixels and consequently the correspondent photoelectrons are added. This process is denominated by binning. The advantages are a faster image acquisition rate, reduced exposition time or an improved signal-to-noise ratio (SNR). However it is achieved at expense of a degradation in image resolution [1]. Spatial sampling rate in object space also depends on objective magnification. For the objectives 40x and 20x these values are 195 nm and 390 nm, respectively, without binning. Consequently putting together e.g., four contiguous pixels, the sampling rate is decreased by the same factor. Spatial sampling in y-axis is imposed by the minimum length in XYZ stage movements in that direction. Based on the specifications of the positioning devices the option was to use the minimum scanning step of 1 μm. On the other hand in microscope applications that demand a larger field-of-view (FOV), wider scanning steps, till the maximum of 20 μm, were used at expense of image lateral resolution.

image n

image i

y-axis

image 1 pixel 1

pixel 1024 2D image x-axis

Figure 6. Schematic of 2D image build. Pixel width is 7,8 μm at image space. Sampling rate in y-axis at object space range from one to twenty microns.

To build a 2D image that is a faithful representation of the relative dimensions of the object in the two lateral directions it is necessary to have an equal scale in both axes. As spatial

A MATLAB-based Microscope http://dx.doi.org/10.5772/58532

sampling rate in the two axes generally is different, the MATLAB function imresize is used. One of its parameters is exactly the ratio of sampling rate values in both axes. This parameter acts as a multiplicative factor to be used for the values in the axis with lower sampling rate. If the 20x objective is used together with the minimum scanning step in y-axis a multiplicative factor of 2,56 should be used in imresize function. Other parameter is the interpolation method that may be chosen from the following three: nearest-neighbor, bilinear or bicubic. This process for 2D image build is illustrated in the scheme of figure 7. It is presented as example one image of the USAF resolution target used in this work for image quality assess‐ ment.

(1)

(2)

(3)

(4)

Figure 7. Diagram depicting the 2D image build in the case of the USAF resolution target. (1) Linear image (data from one sensor readout); (2) 3D representation of the intensity values in the XY plane (after the scanning along y-axis); (3) 2D image using data from sensor readout (no data processing); (4) 2D image with equal scales in both lateral axes (output of imresize function). In (3) and (4): horizontal direction – y-axis; vertical direction – x-axis.

In some microscope applications the field-of-view from each sensor readout is insufficient. However this microscope is able to increase the field-of-view through the scanning also along

57

58

MATLAB Applications for the Practical Engineer

sensor direction (x-axis). The overall width of the 1024 sensor pixels gives a total field-of-view of approximately 400 μm and 200 μm, for 40x and 20x objective lenses, respectively. In order to build the 2D image with a larger field-of-view it is necessary to use an adequate method to mount contiguous images in x-axis. The MATLAB function montage performs this action assuring a correct position alignment. In the overall image these transitions are undistinguish‐ able. 4.2. Extraction and visualization of 3D information Linear images from each sensor readout of its 1024 pixels as well as 2D images obtained in result of the scanning in y-axis are blurred by light collected from reflection on points of the object at different depths. Then this blur in an ideally focused image is a result of light coming from reflection on points in planes in front or behind the focal plane. Scanning optical microscopes (SOM) and particularly the confocal microscope, which is the most widely used, are suitable to get three-dimensional (3D) information. To achieve an image with 3D information the acquisition of different optical sections must be per‐ formed. It consists on images of object planes at different depths. Its acquisition is performed through the scanning along the direction of the optical axis, usually known as the axial direction. Thus as this bench-microscope is intended for the acquisition of images with 3D information, scanning of the XYZ object translation stage is also performed along optical axis (z-axis). Spatial sampling rate in this axis is similar to the used in y-axis. Sampling intervals range from the minimum of 1 μm till 20 μm. The axial resolution depends on the numerical aperture (NA) of the objective. So the selection of the spatial sampling rate must have in consideration which objective is used, namely the 0.4 NA or 0.65 NA. Amongst the different modes of visualization of 3D information the following two are the more usual: • Auto-focus images – images with three-dimensional information as points in image are focused at different depths. • Topographic maps – a three-dimensional representation of the height, h ( xi , yi ), of each point in XY plane. A more detailed description of the methods implemented in this microscope will be presented together with the presentation of several microscope applications. However to illustrate how 3D information can be visualized and simultaneously to evidence the ability of this benchmicroscope to achieve this goal see figure 8. It is the result of the application of a particular non-automatic method. In this case the user select for each lateral position the best focus on the bonding wire amongst the overall axial positions. It is clearly distinguishable the wire inclination that departs from the pad on the integrated circuit and increases its height from the left to the right.

A MATLAB-based Microscope http://dx.doi.org/10.5772/58532

Figure 8. Best-focus image (bonding wire). FOV: 620 μm x 400 μm (horizontal direction – y-axis; vertical direction – xaxis).

5. Results 5.1. Image acquisition The image acquisition as well as its visualization is configured by the computer user. Even in the early stage of this bench-microscope in which preliminary tests were performed with PAF hardware and software it was accomplished. Using PAF system as well as in the initial tests using this sensor readout module the software development tool had been Borland C. This occurred mainly because the software already developed with this tool for PAF was in an easy and fast mode applied to this new hardware module. However a point had been reached in which the need of improving the quality of the graphical user interface and of creating a more user-friendly interface had demanded the change to other software environment. Borland C had been a great tool to develop code for someone that has learned C language in the first step of programming learning. Nevertheless it was not directed for the development of graphical user interfaces (GUI). Consequently it was a time consuming task and the graphical interface had not reached the desired quality. As already mentioned the option had been to use MATLAB. In fact it is a very complete tool as through the use of its objects it is easy and much faster to create a graphical user interface, to add or remove any functionality at any time. But in this bench-microscope there is also the necessity of implement the acquisition of sensor data as well as the control of the positioning of the object stage. Yet MATLAB covers this scope through its toolboxes namely Instrument Control or Data Acquisition. The mode how this functionalities had been easily implemented in MATLAB had been probably the major proof that it was the best choice. CompleteGUI had been the first MATLAB application developed for this bench-microscope. It is the more general in the sense it is used anytime the user intends to get a microscope image. Besides the initialization of the sensor and positioners the user has to define the acquisition parameters in terms of the scanning axes, range and steps.

59

60

MATLAB Applications for the Practical Engineer

As the sensor readout is performed with an acquisition rate of about one frame of 1024 pixels per second and it is a linear image sensor it is necessary to complete the acquisi‐ tion of a set of frames to build a two-dimensional (2D) image. It usually takes a few tens of seconds. This image is then built from a set of one-dimensional (1D) images which are also displayed in real-time. One functionality of this application is to build this 2D image giving the possibility of image visualization immediately after its acquisition is complet‐ ed. This GUI is shown in figure 9 with an example of real-time visualization of raw-data from sensor which is the 1D image.

Figure 9. Graphical User Interface – CompleteGUI – developed for configuration settings and acquisition control and visualization.

The functionalities of CompleteGUI application are summarized below: • Object-stage positioning ◦ User definition of the order of scanning in the three axes (through the GUI); ◦ User definition of the region-of-interest limits in the three axes (through the GUI); ◦ Send commands to positioners through a computer serial port.

A MATLAB-based Microscope http://dx.doi.org/10.5772/58532

• Sensor data readout ◦ Send a command to the sensor readout module to signal acquisition start (through a computer serial port). ◦ Receive the sensor data, 10-bit values of the 1024 pixels (through a computer serial port). • Image visualization ◦ Real-time visualization (preview) of each sensor readout, i.e., a 1D image that shows sensor raw-data; ◦ Visualization of each 2D image as soon as the acquisition of the set of sensor data is completed • Data-files creation ◦ Open Excel data-files ◦ Store data files in the computer hard disk. 5.2. Image quality assessment One of the most important tasks to be performed had been image quality assessment. After the implementation of overall system architecture that includes the optical layout it was necessary to know whether the images obtained with this bench-microscope are in good agreement with theoretical expressions. In fact the performance of a microscope is deter‐ mined by the quality of its output image. There are two different approaches for this assessment: • Quantitative assessment through the determination of the modulation transfer function (MTF) and point spread function (PSF) which are a measure of contrast and resolution, respectively. • Qualitative assessment using images of specific objects with known dimensions [2] Usually only one of these two different methods is selected. It depends on the specificity of the system. Owing to practical reasons a method for quantitative assessment was used. Namely because this is a bright-field microscope in reflection mode and a new method for the deter‐ mination of PSF would be implemented. 5.2.1. Contrast measurement (determination of MTF) Firstly it was measured the MTF which is an usual representation of the performance of imaging systems. It is used as a measure of image contrast. The general definition of MTF is given by the ratio of modulation depth in the output and input of an optical system in function of the spatial frequency when a sinusoidal target is used as input [3]. Modulation depth (M), normally called contrast, has the following definition:

61

62

MATLAB Applications for the Practical Engineer

M=

I max - I min , I max + I min

(1)

where Imax and Imin are the maximum and minimum intensity values, respectively. It is relative either to irradiance emitted by the object or collected in the image, which are system output and input, respectively. Thus MTF definition is:

MTF (x ) =

M output (x )

M sin usoidal (x )

,

(2)

input

whereξ is the spatial frequency (rigorously it is its component along sinusoidal grid direction) From the different methods for MTF determination [4], the choice had been the scanning method. It consists on the measurement of the dependence of the contrast on spatial frequency using a sinusoidal grid as object. Although the MTF concept is applied to sinusoidal grids, for practical reasons, it is much easier to use square grids. A widely used square grid is the USAF (United-States Air Force) resolution target. Then MTF determination was accomplished through contrast measurement for each group / element in images of the USAF target. This target consists on groups of three lines with different densities. The separation of the lines ranges from a minimum of approximately 2 μm (that corresponds to a maximum spatial frequency of 228 lp/mm) to tens of microns (one line pair/mm). In this manner this method consists on the measurement of the system response to an input which is a square and not a sinusoidal grid. This is the definition similar to MTF but designated by contrast transfer function (CTF), in accordance to:

CTF (x f ) =

( ), ( )

M output x f

M square x f input

(3)

whereξ f is the fundamental spatial frequency. To materialize this implementation the MATLAB application USAF_image had been devel‐ oped. Its graphical user interface is presented in figure 10. The functionalities of USAF_image application are summarized below: • Region-of-interest (ROI) ◦ User definition of ROI limits in the two axes (through the GUI); ◦ ROI visualization;

A MATLAB-based Microscope http://dx.doi.org/10.5772/58532

• Image files ◦ Store image files in JPEG format in the computer hard disk; ◦ Parameter acquisition • Selection of a ROI containing one line; ◦ Determination of maximum and minimum value in that ROI; ◦ Determination of line width given in number of pixels using a predetermined threshold value (usually average of maximum and minimum values).

Figure 10. Graphical User Interface – USAF_image – developed for image quality assessment namely contrast meas‐ urements.

CTF values were then calculated for each spatial frequency, using the groups of three lines in USAF target images as shown in figure 11. Different optical configurations had been used in the acquisition of these two images. Images on the left and right correspond to wide-field and line-illumination, respectively. Despite the importance of the illumination mode for the success of this bench-microscope, it is out of the scope of this chapter. Just for completeness the comparison of image contrast on these two illumination modes will be presented.

63

64

MATLAB Applications for the Practical Engineer

(a)

(b)

Figure 11. Images of USAF target showing some elements (sets of three lines) of group 7 (higher density). (a) widefield illumination; (b) line-illumination. Scale bar is 5 μm (horizontal direction – y-axis; vertical direction – x-axis).

This is also the reason why a detailed explanation of how the measured CTF (response to a square input) is converted to MTF (response to a sinusoidal input) is not presented in this chapter. The mathematical expression used for this conversion is:

MTF (x f ) =

CTF é 2n + 1(x f ) ù p N ë û B2 n +1 å 4 n =0 2n + 1

para

xf 2N + 3

ρ0 ,

(8)

where η is a positive scale factor, ρ (q, qobs ) = kq − qobs k is the shortest distance between the mobile robot and the obstacle, and ρ0 is a positive constant that represents the distance of the obstacle’s influence. The graphical representation associated with the obstacle defined by (8) is shown in Figure 3. Hence, the repulsive force, Frep (q), associated with (8), is determined by    q − qo 1 1  − η , for ρ (q, qobs ) ≤ ρ0 , Frep (q) = −∇Urep (q) = (9) ρ (q, qobs ) ρ0 ρ3 (q, qobs )  0, for ρ (q, qobs ) > ρ0 . Finally, it is worth mentioning that the described artificial potential field method depends upon the relative position of the mobile robot to the obstacle, unlike other methods such as that of Krogh [17] in which the potential field is sensitive to the impact time. Furthermore, implementing the presented method requires knowledge of the coordinates of the mobile robot, the obstacle, and the goal, i.e., q = ( x1 , y1 ), qobs = ( xo , yo ), and qm = ( xm , ym ), respectively. In the present work, it is assumed that these coordinates are already known; nevertheless, the employed method can also be extended to sensors that allow the acquisition of these coordinates in real-time.

84

6

MATLAB Applications for the Practical Engineer

Urep(q)

Y

X

Figure 3. Repulsive potential field.

3. Control for obstacle avoidance of WMR This section presents the hierarchical control developed in [18], and analyzed in [26] for a trailer-like vehicle. Using the kinematic model associated with the differentially-driven WMR, an input-output linearization control is proposed that, in conjunction with the artificial potential field method, enables the accomplishment of the obstacle avoidance task. A PI control is then proposed for each DC motor, allowing the WMR to move. Finally, with the ultimate aim of experimentally accomplishing the obstacle avoidance task with a WMR prototype, a hierarchical control is proposed which merges the said controls, similar to the structure presented in [14, 19–23].

3.1. Control of the kinematic model The mobile robot under study is a vehicle comprising two traction wheels, left and right. These two wheels are identical, parallel to each other, non-deformable, and joined by a shaft. The robot also comprises two omnidirectional wheels, front and rear, that ensure the robot platform remains on a plane. Supposing movement is restricted on an XY plane and that there is no wheel slip, existing literature (see [24]) describe the WMR kinematics as given by

( ωr + ω l ) r cos ϕ, 2 ( ωr + ω l ) r y˙ = sin ϕ, 2 ( ωr − ω l ) r , ϕ˙ = 2l x˙ =

(10)

where ( x, y) denotes the position of the mid-point of the shaft that joins the wheels, ϕ is the angle formed by the WMR symmetry axis and the positive X-axis, ωr and ωl are the

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application

85

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application 7 http://dx.doi.org/10.5772/58392 10.5772/58392

angular velocities of the right and left wheels, respectively, r is the wheel ratio, and 2l is the separation between the wheels. Parameters and variables associated with the WMR are shown in Figure 4. This WMR configuration is known as differential traction. l

r

Figure 4. WMR diagram.

Whereas the derivative with respect to time t, associated with the state variables, is denoted by a dot in equations (10) and (12), in the rest of the chapter it is represented explicitly, i.e., d dt . Since the mathematical model described by (10) presents a noninvertible relationship between the controls, (ωr , ωl ), and the outputs, ( x, y), it is not possible to propose a control via input-output linearization. For simplicity, another reference point associated with the WMR that can therefore be considered is the front part which has the coordinates q = ( x1 , y1 ) (see Figure 4). The coordinates of q expressed in terms of x, y, and ϕ are determined by x1 = x + L cos ϕ,

(11)

y1 = y + L sin ϕ,

where L is the distance from the mid-point of the wheel shaft, ( x, y), to the point q in the direction perpendicular to the shaft. Deriving system (11) with respect to time obtains the WMR kinematic model associated with point q = ( x1 , y1 ), which is given by

with A ( ϕ) = 2







x˙ 1 y˙ 1



cos ϕ − L sin ϕ sin ϕ L cos ϕ

= A ( ϕ)

ωr ωl





,

r 2 r 2l

(12)

r 2

− 2lr



.

Since det A ( ϕ) = − Lr 2l 6 = 0, it is clear that an input-output linearization scheme can be

86

8

MATLAB Applications for the Practical Engineer

proposed for (ωr , ωl ) − ( x1 , y1 ). According to [18], an input-output linearization control that allows the WMR to accomplish the obstacle avoidance task and reach the goal can be written as follows: 

ωr ωl



= q

1 Lr f x2 + f y2 + ε υd



L cos ϕ − l sin ϕ L cos ϕ + l sin ϕ

l cos ϕ + L sin ϕ − (l cos ϕ − L sin ϕ)



fx fy



,

(13)

where υd is a desired constant velocity, ε is a constant value close to zero, and f x and f y are the components of Ftotal in directions X and Y, respectively. When there are n obstacles located within the workspace, f x and f y are determined by 

fx fy



=



f x at + f x 1rep + f x 2rep + . . . + f x nrep f y at + f y 1rep + f y 2rep + . . . + f y nrep



,

(14)

with f x at the attractive force component associated with the goal in direction X, and f x 1rep , f x 2rep , . . . , f x nrep the repulsive force components in direction X associated with obstacles 1, 2, . . . , n, respectively. The description of the terms associated with f y is highly similar to that for the terms associated with f x . Here we present the equations explicitly associated with (14) when it is supposed that, within the workspace, there are one and two obstacles, respectively. • One obstacle: In this situation, in accordance with (5) and (7), the attractive potential fields associated with the goal and its force components in directions X and Y, respectively, are given by i 1 h ξ ( x1 − x m )2 + ( y1 − y m )2 , 2 = − ξ ( x1 − x m ),

Uat (q) =

(15)

f x at

(16)

f y at = −ξ (y1 − ym ).

(17)

Whereas, in accordance with (8) and (9), the repulsive potential fields associated with the obstacle and its force components in directions X and Y, respectively, are determined by  h i2 1 1 1 − η , for ρ(q, qobs1 ) ≤ ρ01 , ρ01 ρ(q,qobs1 ) U1rep (q) = 2 0, for ρ(q, qobs1 ) > ρ01 . ih i ( h 1 1 1 η ρ(q,q ) − ρ01 ρ3 (q,q ) ( x1 − xo1 ) , for obs1 obs1 f x 1rep = 0, for ih i ( h η ρ(q,q1 ) − ρ101 ρ3 (q,q1 ) (y1 − yo1 ) , for obs1 obs1 f y 1rep = 0, for

(18) ρ(q, qobs1 ) ≤ ρ01 , ρ(q, qobs1 ) > ρ01 . ρ(q, qobs1 ) ≤ ρ01 , ρ(q, qobs1 ) > ρ01 .

(19)

(20)

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application

87

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application 9 http://dx.doi.org/10.5772/58392 10.5772/58392

Therefore, for a scenario involving one obstacle, (14) is given by 

fx fy



=



f x at + f x 1rep f y at + f y 1rep



.

(21)

• Two obstacles: The attractive force components associated with the goal, f x at and f y at , are determined by (16) and (17), respectively. Whereas the repulsive force components associated with one of the obstacles, f x 1rep and f y 1rep , are determined by (20) and (20), respectively, the second obstacle, which is associated with the repulsive potential field U2rep (q), has the following force components:  h i2 1 η ρ(q,q1 ) − ρ102 , for ρ(q, qobs2 ) ≤ ρ02 , obs2 (22) U2rep (q) = 2 0, for ρ(q, qobs2 ) > ρ02 . ih i ( h η ρ(q,q1 ) − ρ102 ρ3 (q,q1 ) ( x1 − xo2 ) , for ρ(q, qobs2 ) ≤ ρ02 , obs2 obs2 f x 2rep = (23) 0, for ρ(q, qobs2 ) > ρ02 . ih i ( h η ρ(q,q1 ) − ρ102 ρ3 (q,q1 ) (y1 − yo2 ) , for ρ(q, qobs2 ) ≤ ρ02 , obs2 obs2 f y 2rep = (24) 0, for ρ(q, qobs2 ) > ρ02 . Thus, in a scenario two obstacles, (14) adopts the following expression: 

fx fy



=



f x at + f x 1rep + f x 2rep f y at + f y 1rep + f y 2rep



.

(25)

Finally, it is worth mentioning that for each obstacle present within the workspace there will exist a repulsive potential field; as a consequence, the terms associated with f x and f y will be increased depending on the number of obstacles.

3.2. DC motor control In order to execute the obstacle avoidance task experimentally, the angular velocity profiles ωr and ωl , determined by the upper hierarchy (13), must be reproduced by the DC motors associated with the WMR prototype employed in the present work (see [14]). Thus for the right and left angular velocities of the DC motors to approach ωr and ωl , respectively, a PI controller is implemented for each motor. A DC motor mathematical model expressed in terms of the motor shaft speed ̟ is given by di a = u − R a i a − k e ̟, dt d̟ J = −b̟ + k m i a , dt

La

(26)

88

10

MATLAB Applications for the Practical Engineer

where u is the motor armature voltage, i a is the armature current, k e is the back-electromotive force constant, k m is the motor torque constant, L a is the armature inductance, R a is the armature resistance, J is the rotor and load inertia, and b is the viscous friction constant due to both motor and load. Since motor manufacturers do not generally provide all of the parameter values associated with (26), these values can be simply obtained via the reduction of (26) to a first order system that relates ̟ to u. This is accomplished by assuming that L a ≈ 0 in (26). Hence, the simplified model is determined by d̟ = −α̟ + βu. dt

(27)

Characterization of α and β associated with the prototype’s DC motors –two Engel GNM3150s (24 V, 55 W) with G2.6 gearboxes– was previously carried out in [14]. In this latter study it was found that the simplified model (27) of the right and left motors can be respectively expressed by d̟r = −10.20̟r + 5.51ur , dt d̟l = −10.20̟l + 5.99ul , dt

(28)

where ̟r and ̟l are the right and left angular velocity of the motors, and ur and ul represent the armature voltages of the right and left motors, respectively. Hence, a PI control for (27) that achieves ̟ → ̟ ∗ is determined by u = K p e + Ki

Z t 0

edτ,

(29)

with e = ̟ ∗ − ̟,

(30)

where e is the tracking error, ̟ ∗ is the desired angular velocity trajectory, K p is the proportional gain, and Ki is the integral gain. Next, (29) is applied to the DC motors of the WMR.

3.3. Hierarchical control This subsection presents the connection of the control laws developed via a hierarchical control, a block diagram of which is shown in Figure 5. At the upper level, an input-output linearization control was used for the WMR model, generating the desired velocity profiles, ωr and ωl , for the robot wheels to track. These velocity profiles ensure that the WMR moves from the starting point to the goal point while avoiding the obstacles placed in between.

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application

89

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application 11 http://dx.doi.org/10.5772/58392 10.5772/58392

At the lower level, two PI controllers for the DC motors were considered, ensuring that the actual wheel velocities followed the desired velocity profiles generated at the upper level. In accordance with (28), the mathematical models associated with the right and left motors are expressed as follows:

d̟r = −10.20̟r + 5.51ur , dt y r = ̟r ,

(31)

d̟l = −10.20̟l + 5.99ul , dt yl = ̟l .

(32)

and

Thus, two controls, ur and ul , are required; in accordance with (29), these are given by ur = K pr er + Kir ul = K pl el + Kil

Z t 0

Z t 0

er dτ,

(33)

el dτ,

(34)

where ur and ul are the control voltages for the right and left motors, respectively, and K pr , Kir , K pl , and Kil are the constant gains (proportional and integral) associated with each motor. Finally, er and el represent the angular velocity tracking errors defined by

er = ̟r∗ − ̟r , el = ̟l∗ − ̟l ,

(35)

(̟r∗ , ̟l∗ ) = (ωr , ωl ) .

(36)

with
This means that the desired angular velocity trajectories ̟r∗ , ̟l∗ are determined by (ωr , ωl ), which are obtained from (13).

4. Simulations Using Matlab-Simulink, which allows the programming of mathematical models by means of blocks that facilitate the establishment of equations in a transparent and simple manner, this section presents the numerical simulations associated with the obstacle avoidance task developed previously, with the general results then applied to a scenario involving three obstacles within the WMR workspace. As mentioned earlier, the positions of the obstacles

90

12

MATLAB Applications for the Practical Engineer

Figure 5. Block diagram of the WMR hierarchical control.

within the workspace is assumed to be already known. The following subsections outline the stages of the simulation process, namely: the definition of the parameters of the system and of the hierarchical control, the implementation of the control via Matlab-Simulink, and finally the obtained results.

4.1. Definition of parameters When implementing simulations associated with the closed-loop system, both WMR and control parameters must first be considered in order to use them later to program (via Matlab-Simulink) the block associated with the hierarchical control. In this case, the parameters associated with the WMR are l, r, and L, and those associated with the hierarchical control being υd , ε, η, ξ, and ρ0 due to the use of the artificial potential field method. The values of these parameters are shown in Table 1. It is worth mentioning that the values associated with the WMR, i.e., l, r, and L, correspond to the prototype reported in [14]. Finally, the coordinates associated with the goal, qm , and the obstacles, qobs , (for a three obstacle scenario) are presented in Table 2, where ρ0 = ρ01 = ρ02 = ρ03 . The declarations elaborated in Matlab-Simulink for all the above-mentioned parameters are shown in the upper part of Figure 6. Constant l r L υd ε η ξ ρ0

Definition Distance from ( x, y) to the wheels Wheel ratio Distance from ( x, y) to q Desired constant velocity Constant value close to zero Positive scale factor (repulsive) Positive scale factor (attractive) Distance of the obstacle’s influence

Table 1. Parameters employed in the development of the simulations.

Value 0.220 m 0.075 m 0.250 m 0.5 m/s 0.1 2 1 0.5 m

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application

91

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application 13 http://dx.doi.org/10.5772/58392 10.5772/58392

q ( x1 , y1 ) x1 = 0 m y1 = 0 m

qm ( xm , ym ) xm = 2.2 m ym = 1.7 m

qobs 1 ( xo1 , yo1 ) xo1 = 0.2 m yo1 = 0.4 m ρ01 = 0.5 m

qobs 2 ( xo2 , yo2 ) xo2 = 1.2 m yo2 = 0.4 m ρ02 = 0.5 m

qobs 3 ( xo3 , yo3 ) xo3 = 1.5 m yo3 = 1.6 m ρ03 = 0.5 m

Table 2. Coordinates associated with the goal and the obstacles.

Figure 6. Implementation of the hierarchical control in Matlab-Simulink for the WMR.

Once the values of the parameters associated with the WMR and the control have been defined, the hierarchical control can be programmed in Matlab-Simulink to simulate the closed-loop system.

4.2. Implementing the control via Matlab-Simulink Figure 6 (bottom part) shows the three blocks associated with the hierarchical control programmed using Matlab-Simulink: Input-output linearization control, PI controllers/DC motors, and Kinematic model of the WMR. 1.- Input-output linearization control block. In this block, the control determined by (13), which accomplishes the obstacle avoidance task, is programmed. The inputs are the variables ( x1 , y1 , ϕ), and the outputs the desired velocity profiles (ωr , ωl ) = (̟r∗ , ̟l∗ ), the latter

92

14

MATLAB Applications for the Practical Engineer

being tracked by the angular velocities of the right and left DC motors, respectively. Figure 7 presents the block’s constituent sub-blocks. 2.- PI controllers/DC motors block. In this block, the PI controls, associated with the right and left motors, respectively determined by (33) and (34), are implemented. The inputs are determined by (̟r∗ , ̟l∗ ) and the angular velocities produced by the motors (̟r , ̟l ), respectively, with the outputs being the controls ur and ul in such a way that (̟r , ̟l ) → (̟r∗ , ̟l∗ ). The gains associated with the PI controls were here selected as K pr = 2, Kir = 50, K pl = 2, and Kil = 50. 3.- Kinematic model of the WMR block. In this block, the kinematic model of the WMR associated with the point q, determined by (12), is programmed. The inputs are (̟r , ̟l ), and the outputs the variables ( x1 , y1 , ϕ).

Figure 7. Input-output linearization control block.

4.3. Simulation results Finally, this subsection presents the simulation results associated with the obstacle avoidance task for the differentially-driven WMR. The simulations consider the presence of three obstacles within the workspace which have an influence over the WMR during its journey to the goal. Table 2 presents the coordinates associated with the goal and the three obstacles, with the distribution depicted in Figure 8(a) and Figure 14. The results obtained for the variables of interest to the WMR are shown in Figure 8; in this figure one can observe how the WMR successfully evades the obstacles and reaches the goal.

5. Description of the employed prototype This section provides a general description of the WMR prototype reported in [14], which was built to carry out various control tasks associated with WMRs. Figure 9 displays a block diagram describing the existent connections between the different stages of the employed WMR, namely: Subsystems, Power system, and Data acquisition and control system. Here we present a summarized description of the blocks composing the WMR prototype, including its connection with the DS1104 board, as shown in Figure 9. • Stage 1: Subsystems. This stage (see Figure 10) comprises the mechanical subsystems a and b, which include the actuators, sensors, and the mechanical structure of the WMR. Subsystem a, actuators and sensors, generates the movement of the WMR wheels in a specified workspace, and discrete position sensing, respectively. Subsystem b, mechanical

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application

93

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application 15 http://dx.doi.org/10.5772/58392 10.5772/58392

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8. Simulation results in the three obstacles scenario.

94

16

MATLAB Applications for the Practical Engineer

Figure 9. General block diagram of the WMR prototype.

structure, corresponds to the mechanical topology associated with the differentially driven WMR. • Stage 2: Power system. This stage enables interaction between the electronic control interface (stage 3) and the mechanical subsystems (stage 1). A block diagram of stage 2 is shown in Figure 11. This block comprises the following four substages, numbered from 1 to 4: source circuit, optoisolator circuit, H-bridge circuit, and encoder circuit. In substage 1, the power supply, via the source circuit, distributes the different voltages to the general electronic system. Substage 2 enables electrical signal isolation between the DS1104 electronic board and substage 3. Substage 3 enables the direction of rotation of the DC motors to be controlled via the use of a positive or a negative voltage, which is determined by the control exerted by the DS1104 board. Finally, substage 4 involves the acquisition of the encoders’ signals, which can then be used to estimate the position of the WMR within the workspace. • Stage 3: Data acquisition and control system. The main device involved in this stage is the DS1104 board, which performs the acquisition of the variables of interest to the experimental implementation of the WMR hierarchical control. This board was selected due to the potential for integration between Matlab-Simulink and the board’s firmware. Moreover, the high programming level available in Simulink makes it a practical selection for the programming of complex control strategies in a graphical environment. Stage 3 also includes an interface circuit (see Figure 9) which establishes communication between the DS1104 board and the WMR. Figure 12 shows pictures of the real WMR, including the employed instrumentation.

6. Real-time experiments In order to validate the data obtained via numeric simulation presented in Section 4, here we present the experimental results obtained using the WMR in real-time. These experiments were performed using Matlab-Simulink, ControlDesk, and a DS1104 board (dSPACE), with the distribution of the obstacles and the goal identical to that in the simulations.

6.1. Real-time control of the WMR via Matlab-Simulink Whereas the mathematical models associated with the WMR and the DC motors were used to obtain the simulation results, these models were replaced by the WMR prototype in obtaining the experimental results. However, it is worth mentioning that the robot kinematic model

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application

95

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application 17 http://dx.doi.org/10.5772/58392 10.5772/58392

Figure 10. Bottom view of the mechanical structure.

Figure 11. Block diagram of the power system.

was employed here to provide an estimation of the WMR position, since no localization sensor was used for the WMR in the present study. The Matlab-Simulink program associated with the hierarchical control used to obtain the experimental results for the closed-loop system is shown in Figure 13, which consists of the following blocks: System parameters, Input-output linearization control, PI controllers/DC motors, and Kinematic model of the WMR. Figure 13 also depicts the connections between the WMR and the hierarchical control. A comparison of Figures 6 and 13 reveals that the two are highly similar, with the only block experiencing any significant changes being the one associated with the PI controllers/DC motors. Hence, we can focus our attention solely on this block. • PI controllers/DC motors block. This block describes the implementation of the PI controls, ur and ul , associated with the right and left motors, respectively, which allow (̟r , ̟l ) →

96

18

MATLAB Applications for the Practical Engineer

Figure 12. WMR prototype.

(̟r∗ , ̟l∗ ). In order to start the DC motors, the ur and ul voltages must be conditioned through PWM blocks, which are implemented via the DS1104 board. The PWM signals then go through the power system stage illustrated in Figure 9. For the acquisition of the angular velocities ̟r and ̟l , two E50S8 incremental encoders (Autonics) were employed in combination with Matlab-Simulink blocks.

6.2. Experimental results This subsection presents the experimental results associated with the distribution of the obstacles and the goal mentioned in Table 2. Likewise, the locations of the WMR, obstacles, and goal within the workspace are illustrated in Figure 14. The corresponding experimental results for such a scenario, in which the parameters of the WMR and the gains of the controls were the same as those employed in the simulations, are shown in Figure 15. Additionally, supplementary material associated with these results can be seen through Video S1 (see [25]). Also, experimental results for the scenarios involving two and zero obstacles within

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application

97

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application 19 http://dx.doi.org/10.5772/58392 10.5772/58392

Figure 13. Experimental implementation of the controller via Matlab-Simulink and its connection with the WMR.

the workspace are presented in Videos S2 and S3, respectively, (see [25]). The aforementioned videos are available online at: www.controlautomatico.com.mx/trabajos.html

Figure 14. Position and orientation of the starting point of the WMR.

98

20

MATLAB Applications for the Practical Engineer

(a)

(b)

(c)

(d)

(e)

(f)

Figure 15. Experimental results for the three obstacles scenario.

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application

99

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application 21 http://dx.doi.org/10.5772/58392 10.5772/58392

6.3. Discussion of the results Similar results were obtained in the simulation and experimental tests; in both cases the angular velocities of the DC motors, ̟r and ̟l , tracked the desired velocity profiles imposed by the WMR kinematic model, ̟r∗ and ̟l∗ , respectively. Hence, the obstacle avoidance task was successfully carried out by the WMR, as verified experimentally in Figure 15(a). However, it is worth mentioning that an error arose that could only be appreciated visually and which thus does not appear in Figure 15(a): a gap of approximately 10 cm was observed between point q and point qm . This discrepancy likely occurred due to the fact that the position of q within the workspace was calculated indirectly via the WMR kinematic model. It was also observed that the voltages ur and ul , associated with the DC motors, did not surpass the (−24 V, +24 V) voltage interval: this was convenient since the nominal voltage of the employed DC motors is in the range of ±24 V. As a result it can be confirmed that the developed hierarchical control performed successfully.

7. Conclusions Based on the artificial potential field approach, the present work has provided a solution to the obstacle avoidance task for a differentially-driven WMR via the design of a hierarchical control, for which a step-by-step guide has been presented that details the theory, simulation, and experimental implementation. Simulation results were obtained by developing a Matlab-Simulink program, since Simulink provides a graphical environment that facilitates the analysis, design, and construction of dynamic systems. In order to obtain experimental results, a Matlab-Simulink program was again employed, this time alongside ControlDesk and the DS1104 board. A comparison of the simulation and experimental results associated with the obstacle avoidance task revealed the good performance of the developed hierarchical control. However, as mentioned earlier, an error between the WMR and the goal arises that can only be observed visually, likely due to the fact that the position of the robot is calculated via the kinematic model of the WMR. Future work will aim to locate the WMR via a sensor in order to reduce or eliminate said error.

Acknowledgments R. Silva-Ortigoza, H. Taud, M. Marciano-Melchor, and J. A. Álvarez-Cedillo acknowledge financial support from the Secretaría de Investigación y Posgrado del Instituto Politécnico Nacional (SIP-IPN), SNI-México, and the IPN programs EDI and COFAA. The work of C. Márquez-Sánchez was supported by CONACYT and BEIFI scholarships. V. M. Hernández-Guzmán acknowledges financial support from SNI-México. Finally, Rhomina and Joserhamón deserve a special mention from R. Silva-Ortigoza for their moral support and for being the inspiration for his introduction to the thrilling world of robotics.

Author details R. Silva-Ortigoza1 , C. Márquez-Sánchez1 , F. Carrizosa-Corral2 , V. M. Hernández-Guzmán3 , J. R. García-Sánchez1 , H. Taud1 , M. Marciano-Melchor1 , and J. A. Álvarez-Cedillo1

100

22

MATLAB Applications for the Practical Engineer

1 Instituto Politécnico Nacional, CIDETEC, Área de Mecatrónica, Unidad Profesional Adolfo López Mateos, México, DF, Mexico 2 Instituto Tecnológico de Culiacán, Departamento de Metal-Mecánica, Culiacán, Sin, Mexico 3 Universidad Autónoma de Querétaro, Facultad de Ingeniería, Querétaro, Qro, Mexico

References [1] R. Silva-Ortigoza, M. Marcelino-Aranda, G. Silva-Ortigoza, V. M. Hernández-Guzmán, M. A. Molina-Vilchis, G. Saldaña-González, J. C. Herrera-Lozada, and M. Olguín-Carbajal, “Wheeled mobile robots: A review,” IEEE Latin America Transactions, vol. 10, no. 6, pp. 2209–2217, 2012. Available at http://www.ewh.ieee.org/reg/9/etrans/eng/ [2] J. L. Crowley, “World modeling and position estimation for a mobile robot using ultrasonic rangin,” in Proceedings of the IEEE International Conference on Robotics and Automation, pp. 674–680, May 1989. [3] J. T. Schwartz and M. Sharir, “On the piano movers problem. II. General techniques for computing topological properties of real algebraic manifolds,” Advances in Applied Mathematics, vol. 4, no. 3, pp. 298–351, 1983. [4] N. J. Nilsson, “A mobile automaton: An application of artificial intelligence techniques,” in Proceedings of the International Joint Conference on Artificial Intelligence, pp. 509–520, May 1969. [5] O. Khatib, “Real-Time obstacle avoidance for manipulators and mobile robots, ” The international Journal of Robotics Research, vol. 5, no. 2, pp. 90–98, 1986. [6] J. Borenstein and Y. Koren, “Obstacle avoidance with ultrasonic sensors,” IEEE Journal of Robotics and Automation, vol 4, no. 2, pp. 213–218, 1988. [7] J. Borenstein and Y. Koren, “Real-Time obstacle avoidance for fast mobile robots, ” IEEE Transactions on Systems, Man and Cybernetics, vol. 19, no. 5, pp. 1179–1187, 1989. [8] R. Tilove, “Local obstacle avoidance for mobile robots based on the method of artificial potentials,” in Proceedings of the IEEE International Conference on Robotics and Automation, pp. 566–571, May 1990. [9] J. Borenstein and Y. Koren, “Potential field methods and their inherent limitations for mobile robot navigation,” in Proceedings of the IEEE International Conference on Robotics and Automation, pp. 1398–1404, April 1991. [10] J. Guldner and V. I. Utkin, “Sliding mode control for gradient tracking and robot navigation usig artificial potential fields,” IEEE Transactions on Robotics and Automation, vol. 11, no. 2, pp. 247–254, 1995.

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application

101

Obstacle Avoidance Task for a Wheeled Mobile Robot – A Matlab-Simulink-Based Didactic Application 23 http://dx.doi.org/10.5772/58392 10.5772/58392

[11] K. Valanavis, T. Hebert, R. Kolluru, and N. Tsourveloudis, “Mobile robot navigation in 2-D dynamic enviroments using an electrostatic potential field,” IEEE Transactions on Systems, Man and Cybernetics, Part A: Systems and Humans, vol. 30, no. 2, pp. 187–196, 2000. [12] F. Xu, H. Van Brussel, M. Nuttin, and R. Moreas, “Concepts for dynamic obstacle avoidance and their extended application in underground navigation,” Robotics and Autonomous Systems, vol. 42, no. 1, pp. 1–15, 2003. [13] W. Huang, B. Fajen, J. Fink, and W. Warren, “Visual navigation and obstacle avoidance using a steering potential function,” Robotics and Autonomous Systems, vol. 54, no. 4, pp. 288–299, 2006. [14] R. Silva-Ortigoza, C. Márquez-Sánchez, M. Marcelino-Aranda, M. Marciano-Melchor, G. Silva-Ortigoza, R. Bautista-Quintero, E. R. Ramos-Silvestre, J. C. Rivera-Díaz, and D. Muñoz-Carrillo, “Construction of a WMR for trajectory tracking control: Experimental results,” The Scientific World Journal, vol. 2013, Article ID 723645, pp. 1–17, 2013. Available at http://dx.doi.org/10.1155/2013/723645 [15] V. M. Hernández-Guzmán, R. Silva-Ortigoza y R. V. Carrillo-Serrano, Control Automático: Teoría de Diseño, Construcción de Prototipos, Modelado, Identificación y Pruebas Experimentales, Colección CIDETEC–IPN, Mexico, DF, Mexico, 2013. Available at http://www.controlautomatico.com.mx [16] S. S. Ge and Y. J. Cui, “New potential functions for mobile robot path planning,” IEEE Transactions on Robotics and Automation, vol. 16, no. 5, pp. 615–620, 2000. [17] B. H. Krogh, “A generalized potential field approach to obstacle avoidance control,” International Robotics Research Conference, Bethlehem, PA, August 1984. [18] J. R. García Sánchez, Diseño y construcción de un robot móvil, aplicando el método de campos potenciales en la evasión de obstáculos. Tesis de Maestría. CIDETEC del Instituto Politécnico Nacional, Mexico City, Mexico, 2008. [19] A. W. Divelbiss and J. T. Wen, “Trajectory tracking control of a car-trailer system,” IEEE Transactions on Control Systems Technology, vol. 5, no. 3, pp. 269–278, 1997. [20] R. Silva-Ortigoza, G. Silva-Ortigoza, V. M. Hernández-Guzmán, V. R. Barrientos-Sotelo, J. M. Albarrán-Jiménez, and V. M. Silva-García, “Trajectory tracking in a mobile robot without using velocity measurements for control of wheels,” IEEE Latin America Transactions, vol. 6, no. 7, pp. 598–607, 2008. Available at http://www.ewh.ieee.org/reg/9/etrans/eng/ [21] R. Silva-Ortigoza, J. R. García-Sánchez, J. M. Alba-Martínez, V. M. Hernández-Guzmán, M. Marcelino-Aranda, H. Taud, and R. Bautista-Quintero, “Two-stage control design of a Buck converter/DC motor system without velocity measurements via a Σ − ∆-modulator,” Mathematical Problems in Engineering, vol. 2013, Article ID 929316, pp. 1–11, 2013. Available at http://dx.doi.org/10.1155/2013/929316

102

24

MATLAB Applications for the Practical Engineer

[22] R. Silva-Ortigoza, C. Márquez-Sánchez, F. Carrizosa-Corral, M. Antonio-Cruz, J. M. Alba-Martínez, and G. Saldaña-González, “Hierarchical velocity control based on differential flatness for a DC/DC Buck converter-DC motor System,” Mathematical Problems in Engineering, vol. 2014, Article ID 912815, pp. 1–12, 2014. Available at http://dx.doi.org/10.1155/2014/912815 [23] R. Silva-Ortigoza, V. M. Hernández-Guzmán, M. Antonio-Cruz, and D. Muñoz-Carrillo, “DC/DC Buck power converter as a smooth starter for a DC motor based on a hierarchical control,” IEEE Transactions on Power Electronics, Article in Press, Available with DOI: 10.1109/TPEL.2014.2311821 [24] H. Sira-Ramirez and S. K. Agrawal, Differentially Flat Systems, Marcel Dekker, New York, 2004. [25] Supplementary Material. Available at http://www.controlautomatico.com.mx/trabajos.html [26] T. A. Vidal Calleja, Generalización del método de campos potenciales artificiales para un vehículo articulado. Tesis de Maestría. Sección de Mecatrónica del Departamento de Ingeniería Eléctrica del CINVESTAV-IPN, Mexico City, Mexico, 2002.

Chapter 4

Dual Heuristic Neural Programming Controller for Synchronous Generator Mato Miskovic, Ivan Miskovic and Marija Mirosevic Additional information is available at the end of the chapter http://dx.doi.org/10.5772/58377

1. Introduction This chapter describes the development of voltage control system of a synchronous generator based on neural networks. Recurrent (dynamic) neural networks (RNN) are used, as a type that has great capabilities in approximation of dynamic systems [1]. Two algorithms are used for training – Dual Heuristic Programming (DHP) and Globalized Dual Heuristic Program‐ ming (GDHP). The algorithms have been developed for the optimal control of nonlinear systems using dynamic programming principles. Neural voltage controller is developed in MATLAB and Simulink environment. For training purposes a mathematical model of synchronous generator is designed and applied in Simu‐ link. DHP and GDHP algorithms are designed in Simulink, with matrix calculations in Sfunctions. Algorithms are used for offline training of neural networks (NN). In the second part, the same functions are redesigned as real time controllers, based on the Real Time Windows Target Toolbox. DHP or GDHP algorithms provide a significant improvement over conventional PI controller with, in some cases, power system stabilizer (PSS). Conventional linear controller can only provide optimal control in single operating point. Algorithms described in the chapter provide significantly better control over the whole operating range by minimization of the user defined cost function during the training of the neural network. Digital voltage controller is implemented on real system in two basic steps. First step is neural network design and training in Matlab environment on desktop computer. Second step consists of transfer of controller with computed coefficients to the digital control system hardware.

104

MATLAB Applications for the Practical Engineer

For the controller to have optimal performance in all real operating conditions, neural networks must be trained for the whole working range of the plant. Procedure described in the chapter can be applied on standard voltage control systems by replacing the PI controller with the offline-trained neural network. Most modern digital control systems have sufficient hardware and software support for neural network implementation.

2. Recurrent neural networks Neural networks are computational models capable of machine learning and pattern recog‐ nition. They can also be used as universal approximators of continuous real functions. Typical dynamic neural network is shown in Figure 1.

Σ Σ

u1 (t )

ϕ

w o ij

w hi j

w h i

j

u2 (t )

Σ Σ

ϕ

Σ Σ

ϕ

x(t)

Σ Σ

Y1 (t )

ϕ

x (t − 1) wij

Figure 1. Recurrent neural network structure

Figure 1. Recurrent neural network structure

Neural networkfrom is formed from interconnected layers layers of neurons. Output of Output each neuronof is each neu network is formed interconnected of neurons. calculated from: ted from:  m    æ m ö  y =y ϕ= j ç å∑( w xw) + bx÷  + b  (1) k k çè j = 0 kj kj kj ÷ø kj   = 0 j  

ϕ is activation function,

w

kj

and x

kj

are weights and outputs from previous layer and b is bias.

activation function maps input values

u (t ) ∈ −∞, ∞ 

to output

y (t ) ∈ 0,1

or

y (t ) ∈

Dual Heuristic Neural Programming Controller for Synchronous Generator http://dx.doi.org/10.5772/58377

where

φ is activation function, wkj and xkj are weights and outputs from previous layer and b is bias. Typical activation function maps input values u(t) ∈ − ∞, ∞ to output y(t) ∈ 0, 1 or y(t) ∈ − 1, 1 . Activation function must be derivable over the whole range. Among typically used nonlinear activation functions are logsig, tansig and Gaussian. Dynamic behavior is added by introducing a memory element (step delay) that stores neuron output and reintroduces it as input to all neurons in the same layer in next time step. This type of network is named recurrent or dynamic neural network. An example of tansig 2 − 1, and its (Hyperbolic Tangent Sigmoid Transfer Function) activation function, ϕ (u ) = 1 + e −2u derivative ϕ (u )′ = 1 − ϕ (u )2 is shown in Figure 2.

1

Hyperbolic tangent sigmoid transfer function

 (u ) =

2 1 1 + e 2u

0.5

 (u )

0 -0.5 -1 -5

-4

-3

-2

-1

0

1

2

3

4

5

1

2

3

4

5

u 1

Hyperbolic tangent sigmoid transfer derivative function

 (u ) ' =

4e 2u (1 + e 2u ) 2

 (u ) ' = 1   (u )

2

0.8

 ' (u )

0.6 0.4 0.2 0 -5

-4

-3

-2

-1

0

u

Figure 2. Tansig activation function and its derivative

Input values are in range of u ∈ − ∞, ∞ , with output values in range φ(u) ∈ − 1, 1 . Use of function with simple derivative significantly improves calculation times of numerical procedures. Dynamic (recurrent) neural networks provide good approximation of time-dependent functions. However, training of recurrent neural networks is more complex than training static networks.

105

106

MATLAB Applications for the Practical Engineer

2.1. Recurrent neural networks training There are many methods for training neural networks, and most of them rely on some form of backpropagation – calculation of partial derivative of system output over network weight ∂y coefficients . Partial derivatives are determined using chain rule, by finding derivatives ∂w of weight coefficients in output layer, using intermediate results to calculate next layer, with respect to network structure. Calculation is repeated until all layers are processed.

( )

Described procedure is valid for static neural networks, but is of no use in recurrent networks, as each neuron has inputs from previous time step from the same layer. Two procedures are often used in calculation of partial derivatives: Backpropagation Through Time (BPTT) and Real Time Recurrent Learning (RTRL). In BPTT recurrent neuron inputs are replaced with the same neural network delayed by one time step. Process is iterated N times. Thus, recurrent network is approximated by N feedfor‐ ward networks with delayed outputs. Network weight coefficients, recurrent inputs and ∂ y (t ) outputs must be preserved for past N time steps. Obtained derivative is used to ∂ w (t ) determine new weight coefficients. One advantage of BPTT is that all methods used for training feedforward networks can be utilized in training recurrent networks.

(

)

This chapter will describe use of RTRL procedure in training dynamic networks. Output of recurrent network with one hidden and one output layer has the following form:

(

x (t ) = tanh W1 × p1 y (t ) = W2 × x (t )

)

(2)

where

p1– hidden layer input, p1 = u(t) ; x(t − 1) ; 1 , x(t)– hidden layer output, y(t) – network output,

W 1– hidden layer weight coefficients, W 2– output layer weight coefficients, p2 = x(t) ; 1 Equation (2) can easily be expanded for multi layer recurrent networks. Partial derivative of network output over output layer weight coefficients is ¶y (t ) = p2 (t ) ¶wijo (t )

(3)

Dual Heuristic Neural Programming Controller for Synchronous Generator http://dx.doi.org/10.5772/58377

For hidden layer, partial derivatives are ¶y (t ) ¶y (t ) ¶x (t ) = ¶wh (t ) ¶w2 (t ) ¶w2 (t )

(4)

Partial derivative of hidden layer output over weight coefficients is

NW 1 æ ¶xl (t ) ¶xk (t  1) ö = f [ x (t )] 'i ç p j (t )d li + å w(t )l , N in + k ÷ ç ¶wi , j (t ) ¶wi , j (t  1) ÷ø k =1 è

(5)

where

N W 1– number of neurons in hidden layer, N in – number of neurons in input layer. Value of

∂ xl (t) is calculated using four nested loops (Appendix A). ∂ wi, j

Matrix

of

output

derivatives

over

weight

coefficients

is

of

the

form

N W 1, N W 1 ⋅ ( N in + N W 1 + 1) . 2.2. Neural network training algorithms 2.2.1. Gradient descent Gradient descent is the most commonly used method in neural networks training. Difference between real and desired network output can be expressed as e(t ) = y (t )  yˆ (t )

(6)

Differentiation of (5) provides update values of weight coefficients:

wij (t ) = wij (t  1)  b

¶yˆ(t ) ¶wij (t )

Parameter β defines learning rate and must be in range 0 < β < 1.

(7)

107

108

MATLAB Applications for the Practical Engineer

2.2.2. Kalman filter based training Kalman filter [5] was developed as a recursive procedure of optimal filtering of linear dynamic state-space systems. For RNN training extended Kalman filter (EKF) is used. It is also appli‐ cable to nonlinear systems using linearization around operating point. The algorithm is used as a predictive-corrective procedure. Kalman filter based recurrent network training is described in [6] and [7]. Training is performed in two steps. The first step is calculation of prediction values of weight coefficients. Second step is correction of predicted variables based on measured state space. Algorithm is defined by the following relations:

(

K (t ) = P (t ) H (t ) h (t ) I + H T (t ) P (t ) H (t )

)

(8)

W (t ) = W (t  1) + K (t )e (t )

(9)

P (t + 1) = P (t )  K (t ) H T (t ) P (t ) + Q (t )

(10)

where

K (t)– Kalman gain matrix, determined form network output derivatives over weight coeffi‐ ∂ y(t) . System linearization is accomplished by calculating H (t). cients H (t) = ∂w P(t)– Error covariance matrix, calculated in every iteration To escape local minima during training, a diagonal matrix Q(t) is added to covariance matrix P(t) in each time step. Updates to weight coefficients are determined from (9), as a product of Kalman gain and ε(t) in current step. Training based on Kalman filter is fast, has good rate of convergence and is often the only practical solution for recurrent networks training. However, the procedure is demanding in terms of computational power and memory requirements.

3. System identification Nonlinear dynamic systems are identified as input-output system model, or as a state space system. For the selected process, identification is performed by minimization of the cost function

Dual Heuristic Neural Programming Controller for Synchronous Generator http://dx.doi.org/10.5772/58377

E (t ) =

1 T × [ y (t )  yˆ (t )] × [ y (t )  yˆ (t )] 2

(11)

where y(t) − ^y (t) is the difference of the measured system outputs and output values of the model. Selected cost function represents squared error of the identified model. Dynamic system identification can be performed according to [9], using either series-parallel or parallel structure.

u (t )

Non linear dynamic system

y (t )

u (t )

Non linear dynamic system

Er Er

ro

r

ro

r

y (t )

 -

 -

z

z 1 1

Neural Networks

yˆ(t )

Neural Networks

yˆ(t ) z 1

z 1 Series- parallel Identification

Parallel Identification

Figure 3. Nonlinear dynamic system identification

In series-parallel structure Figure 3, model neural network input vector is expanded with system outputs from previous steps. In parallel structure Figure 3, network output from previous states is added to input vector. Series-parallel identification provides estimated system output for next time step only, while parallel identification can predict multiple future iterations. However, parallel structure is harder to train and is susceptible to state drifting. To verify the validity of the RTRL algorithm, neural networks were trained on several discrete mathematical functions. Supervised training has been used, with control value expressed as difference between desired and actual network output. A

parallel identification structure was used on dynamic SISO system y(t) y(t + 1) = + u(t)3. Training has been performed on Recurrent Neural Network with 1 + y(t)2 four neurons in hidden layer, and one neuron in linear output layer. Identification results are shown in Figure 4.

109

MATLAB Applications for the Practical Engineer

(a)

(c) 300

40

250

20

200 E(t)

u(t) y(t)

60

0 -20

-60 0

20

40

60

80

100 k

120

140

160

180

150 100

u(t) y(u) y(u) RNN

-40

50

200

0

0

0.5

1 k

1.5

2 4

x 10

(b)

100

u(t) y(t) y(t) RNN

50

u(t) y(t)

110

0

-50

-100

0

20

40

60

80

100 k

120

140

160

180

200

Figure 4. Identification of dynamic function 1

Figure 4a compares initial response of model output y(t) and network output ^y (t). Figure 4b represents the same values with trained network. There is a perfect match between function output values and output of the trained NN. Figure 4c displays training rate as error square ( y(t) − ^y (t))2. Training was completed in 20000 steps. Input signal used for training and verification is u(t) = sin(

(T D = 1s ).

2π 2π ) + sin( ), with sample time 50 12

In second test, a system with one input, two state variables and one output has been identified: x1(t + 1) = 0.5x2 ⋅ (t) + 0.2 ⋅ x1(t) ⋅ x2(t) x2(t + 1) = − 0.3 ⋅ x1(t) + 0.8 ⋅ x2(t) + u(t) y(t + 1) = ( x1(t) + x2(t))2 System was identified using series-parallel structure with RNN with 5 neurons in hidden layer, and one linear neuron in output layer. Identification results are shown in Figure 5.

Dual Heuristic Neural Programming Controller for Synchronous Generator http://dx.doi.org/10.5772/58377

(a)

120

(c) 350 u(t) y(t) y(t) RNN

100

300 250

60

E(t)

u(t) y(t)

80

40

200 150 100

20

50

0 0

50

100

150

200

250

0

0

0.5

k

1.5

2 4

x 10

(b)

120

u(t) y(t) y(t) RNN

100 80 u(t) k(t)

1 k

60 40 20 0 0

50

100

150

200

250

k

Figure 5. Identification of dynamic function 2

Input signal used for training and verification is the same as in previous example: u(t) = sin(

2π 2π k T D ) + sin( kT D ), where (T D = 1s ). 50 22

Both examples were designed in Matlab and Simulink. Neural network was trained using clanguage S-function (Appendix B). Although the training took 10000 steps, all important system dynamics were visible on trained network after 2000 steps. 3.1. Determination of partial derivatives of functions using neural networks One of the most important features of neural network modeling is the possibility to determine ∂ yi (t) . After finding output derivatives, various system output derivatives over inputs ∂ uj (t)

(

)

optimization algorithms can be used. Partial derivative is calculated from (2) p ¶yi (t ) NW 1 æ = å ç Wi , j × å f k '( x ) × Wk , j ç ¶u j (t ) j =1 è k =1

N

(

ö

) ÷÷ ø

(12)

111

MATLAB Applications for the Practical Engineer

)

(

where N p = N W 1 + N in + 1 is number of hidden layer inputs. Partial derivative

∂ yi (t) is ∂ uj (t)

obtained from network backpropagation.

∂ y(t) was first calculated ∂ u(t) analytically, and then obtained from model network using (12). Results are compared in Figure 6a and 6b.

For nonlinear function y(t) = u(t − 1)2 + 6 ⋅ u(t − 1)3, partial derivative

(a) 3.5 3

u(k) Y(k) y(k) NN

2.5 2 u(k) Y(k)

1.5 1 0.5 0 -0.5 -1

3

3.01

3.02

3.03

3.04

3.05 k

3.06

3.07

3.08

3.09

3.1 4

x 10

(b) 9 dY/du dY/du RNN

8 7 6 5 dY/du

112

4 3 2 1 0 -1

3

3.01

3.02

3.03

3.04

3.05 k

3.06

3.07

3.08

3.09

3.1 4

x 10

Figure 6. Identification of nonlinear function and calculation of partial derivatives using neural networks

Results from Figure 6b show good match between calculated function derivatives and values obtained from neural network of identified model. In S-function for training RNN (Appendix B) a procedure is shown for online calculation of ∂ yi (t) partial derivatives of , using (12). ∂ uj (t)

(

)

4. Adaptive critic design Adaptive critic design (ACD) is a set of optimal control procedures of nonlinear systems. It is based on dynamic programming and neural networks. Optimal control is achieved by training

Dual Heuristic Neural Programming Controller for Synchronous Generator http://dx.doi.org/10.5772/58377

neural networks and using them in the structure according to the selected ACD algorithm. ACD algorithms are implemented using a heuristic approach – system state space variables, control values and cost function are selected based on process characteristics. This text focuses on the Heuristic Dynamic Programming (HDP), Dual Heuristic Programming (DHP) and General Dual Heuristic Programming (GDHP) [1,2]. 4.1. Heuristic dynamic programming Heuristic dynamic programming (HDP) is the simplest (by structure) ACD algorithm. The basis of all HDP algorithms is minimization of Bellman dynamic programming function J (t): ¥

J (t ) = å g kU (t + k ) k =0

(13)

where

γ represents discount factor with value in range 0 < γ < 1, ensuring the convergence of Bellman sum J (t) and cancellation of old values, U (t + 1) is user defined utility function based on system state space variables. Control signal that minimizes Bellman sum is generated using one of the Adaptive Critic algorithms. In practical control systems, HDP optimization is realized with three neural networks with a structure as shown in Figure 7. Model neural network in 7a is used for model identification. Network inputs are control signal u(t) and output state in the previous time step y(t − 1). Network output represents estimated model output in the current time step ^y (t). Only systems outputs used in utility function need to be estimated. Model NN is structured according to the selected process model, where the control value is input to the NN, and NN output represents predicted system output. The same NN output values are used to calculate the cost function. Model neural network is trained to minimize the quadratic criterion: E M = å E TM (t ) × E M (t )

(14)

E M = y (t )  yˆ (t )

(15)

t

where:

For adaptive control systems, model NN is trained continuously, allowing for adaptation to change of system parameters, and optimal control. Model neural network is trained simulta‐ neously with other networks over time. Critic neural network in 7b is used for identification of Bellman sum J (t).

113

114

convergence of Bellman sum J (t ) and cancellation of old values, U (t + 1) is user defined utility function based on system state space variables. Control signal that minimizes Bellman sum is generated using one of the Adaptive Critic algorithms. In practical control systems, HDP optimization is realized with three neural MATLAB Applications the Practical Engineerin Figure 7. networks with aforstructure as shown

yˆ(t + 1)

u (t ) TD

y (t − 1)

TD

J (t + 1)

TD

yˆ(t )

MODEL NETWORK

CRITIC TD

TD TD

J ( t ) − [γ ⋅ J (t + 1) + U (t ) ]

yˆ(t ) y ( t ) − yˆ ( t )

y (t )

TD

CRITIC

J (t )

TD

(a)

(b)

u (t ) y (t )

TD

y (t − 1)

y(t) TD

MODEL NETWORK

uˆ(t )

yˆ(t ) TD

ACTION TD

TD

∂y (t ) ∂u(t )

TD

X

∂J (t ) ∂y (t ) TD

CRITIC TD

J (t )

(c)

Figure 7. HDP Control structure

Figure 7. HDP Control structure

Model neural network in 7a is used for model identification. Network inputs are control signal

Training is used tostate minimize output in the previous time step y (t − 1) . Network output represents estimated u (t ) and model output in the current time step yˆ(t ) . Only systems outputs used in utility function need EC = å ECT (t ) × EC (t )

(16)

EC (t ) = J (t )  [g × J (t + 1) + U (t )]

(17)

t

where:

Critic neural network training is performed with a copy of the same network which is fed with predicted future system outputs ^y (t + 1), thus approximating J (t + 1). Action neural network training is based on model and critic network outputs, by minimizing criterion: E A = å E TA (t ) × E A (t ) t

(18)

Dual Heuristic Neural Programming Controller for Synchronous Generator http://dx.doi.org/10.5772/58377

where: E A (t ) =

Partial derivatives

¶y (t ) ¶J (t ) × ¶u (t ) ¶y (t )

(19)

∂ y(t) ∂ J (t) and are determined from model network using backpropa‐ ∂ u(t) ∂ y(t)

gation in the manner described in section 2.5. Minimization of scalar product

EA in (18) by training of action network results in minimi‐

zation of Bellman sum in J (t) in (13). Minimum of J (t) is achieved when

∂ J (t) ≈ 0, as action ∂ y(t)

network training criterion approaches zero, thus ensuring end of training process. Action neural network is trained from the output network (unsupervised algorithms). Training of neural networks is carried out using (19) with the training concept shown in Figure 7c. 4.2. Dual heuristic programming Dual Heuristic Programming optimal control algorithm is developed with a goal of increasing procedure efficiency by increasing training speed. Model network training is the same as in HDP procedure. The improvement is in the direct calculation of the Bellman sum partial derivative over state space variables

∂ J (t + 1) as output ∂ y(t)

of critic network. Structure of DHP control system is shown in Figure 8. Training of critic neural network minimizes scalar product: EC = å ECT (t ) × EC (t ) t

(20)

with

EC (t ) =

where

¶J [ yˆ (t )] ¶J [ yˆ (t + 1)] ¶U [ y (t )] g  ¶yˆ(t ) ¶y (t ) ¶y (t )

(21)

∂ J ^y (t) represents partial derivative of Bellman sum over state space variables, and ∂ ^y (t)

is estimated from critic network. The second part of (21) is determined by deriving

115

116

procedure efficiency by increasing training speed. Model network training is the same as in HDP procedure. The improvement is in the direct ∂J (t + 1) calculation of the Bellman sum partial derivative over state space variables as ∂y (t ) MATLAB Applications for the Practical Engineer output of critic network. Structure of DHP control system is shown in Figure 8. yˆ(t + 1)

u(t ) TD

y (t − 1)

TD

MODEL NETWORK

yˆ( t )

TD

TD

λ (t + 1) ⋅

TD

λ (t + 1) CRITIC

TD

∂y (t + 1) ∂y ( t ) + -

∂y (t + 1) ∂u(t + 1) λ (t + 1) ⋅ ⋅ ∂y (t ) ∂y (t )

-

∂U (t ) ∂u (t )

∂U (t ) ∂u(t ) ⋅ ∂u (t ) ∂y (t )

y (t )

yˆ( t )

ACTION

TD TD

-

CRITIC

TD

λ (t )

TD

uˆ(t ) (a)

yˆ(t + 1)

u(t ) TD

y (t − 1)

TD

MODEL NETWORK

λ (t + 1) CRITIC

TD

yˆ( t )

TD

TD TD

λ (t + 1) ⋅

∂y (t + 1) ⋅γ ∂u (t ) +

∂U (t ) ∂u (t )

+

y (t ) ACTION

TD TD

(b)

Figure 8. DHP structure, (a) Critic NN training, (b) action NN training

Figure 8. DHP structure, (a) Critic NN training, (b) action NN training

¶J ( t + 1) n ¶yˆ ( t + 1) m n ¶yˆ ( t + 1) ¶uˆk ( t + 1) = å li ( t + 1) i + åå li ( t + 1) i ¶yˆ j ( t ) ¶yˆ j ( t ) ¶uˆk ( t ) ¶yˆ j ( t ) i =1 k =1 i =1

where ^ ^ (t + 1) = ∂ J y (t + 1) , λ ∂ ^y (t + 1) n is size of model network output vector ^y (t),

^ (t). m is size of control vector u

(22)

Dual Heuristic Neural Programming Controller for Synchronous Generator http://dx.doi.org/10.5772/58377

Third part of equation (22) is determined from deriving ¶U ( t ) ¶Y ( t )

m

=å

¶U ( t ) ¶uk ( t + 1)

(23)

k =1 ¶uk ( t ) ¶y j ( t )

∂ ^y i (t + 1) , ∂ ^y j (t ) ^ ^ ∂ y i (t + 1) ∂ U (t ) ∂ u k (t + 1) are determined from model network, and partial derivatives , ^ ( ) ( ) ∂ uk t ∂ ^y j (t ) ∂ uk t ^ from action network. Prediction of states y (t + 1) and λ(t + 1) is determined from model network. Somewhat simplified structure of DHP is shown in Figure 8a. Partial derivatives

Action neural network is trained to minimize Bellman sum, min J (t) , from the criterion ∂ J ^y (t) ≈ 0. ∂ ^y (t) Thus, equation (22) can be transformed to g

¶J [ yˆ(t + 1)] ¶y (t )

+

¶U [ y (t )] ¶y (t )

(24)

=0

Training of Action Network minimizes scalar product: E A = å E TA (t ) × E A (t )

(25)

t

with E A (t ) = g

¶J [ yˆ(t + 1)] ¶y (t )

+

¶U [ y (t )]

(26)

¶y (t )

Simplified training procedure of action NN is shown in Figure 8a. Cross training of critic and action networks results in slower approximation. In order to reduce training time, it is recommended to begin the procedure with action network pre-trained as linear controller. 4.3. Global dual heuristic programming Globalized Dual heuristic optimal control algorithm uses both algorithms described. Critic neural network provides approximation of Bellman sum J (t) and derivative

∂ J (t + 1) ∂ y(t)

117

118

MATLAB Applications for the Practical Engineer

Action network is trained using the same procedure as described for HDP and DHP.

5. Neural networks based control of excitation voltage Main elements of electric power system are generators, transmission and distribution network and consumers. Transmission and distribution network form a power grid, which is usually very large in geographic terms. Grid control is implemented on both local and global level, mostly on the power generation side of the grid. Specific to power grid control is large number of closed control loops governing the power and voltage of synchronous generators. Trans‐ mission network is governed by setting voltage on transformers in order to minimize trans‐ mission losses. Governing of generator power and voltage achieves balance of power as well as reactive power flow on the transmission system. Voltage control also has a large impact on the dynamic system behavior. Generator power output is governed on the generator unit. For water and steam turbine generators, it defines local power transmission to the grid, while also influencing system frequency. Considering the inertness of the generator systems, power output is not difficult to govern. Voltage control manages the balance between generator and grid voltage, and determines reactive power flow. A change in generator voltage setpoint has impact on system state variables and parameters. Disregarding system transients, a change of generator voltage causes a change on generator busbars, generator current and reactive power. If we also consider system dynamics, generator voltage change also changes system parame‐ ters influencing system damping. Synchronous generator connected to the power grid is dominantly nonlinear in characteristics. Generator output power is proportional to generator and grid voltage, while the relation to angle between generator and grid voltage is a sinus π function. Maximum theoretical output power is realized with a load angle being δ = . 2 Exceeding of maximum value of load angle causes large disturbances and disconnection from power grid. Operating point of generator connected to the grid is determined by steady state excitation current. Excitation voltage and current values dictate generator voltage, reactive power and load angle. Frequent small disturbances expected and part of normal operating conditions, resulting in oscillations of generator electromechanical variables. Especially pronounced is the swinging of generator power, reaching several percents of nominal value, with a frequency around 1Hz. It is caused by low system damping and outer disturbances on the grid. Also present on the grid are oscillations caused by the system characteristic frequencies, usually around 0.1Hz. Power grid can be represented as a single machine connected to an infinite bus system. For each operating point a characteristic frequency and damping can be determined. In case of small disturbances with a synchronous generator working around predefined operating point, system can be considered linear, but only if it returns to the same operating

Dual Heuristic Neural Programming Controller for Synchronous Generator http://dx.doi.org/10.5772/58377

point after each disturbance. For those conditions, power system stabilizer (PSS) is added to AVR. Based on change of output power and load angle, it modifies voltage regulator setpoint and increases system damping. PSS is configured using linear control theory, and is presently a standard part of automatic voltage regulator systems [10]. Large disturbances occur after short circuits on the power grid, transmission and generator outages, and connection and disconnection of large loads. Electromechanical oscillations occurring as a consequence of disturbances cause damped oscillations of all electrical and mechanical values (generator voltages and currents, active and reactive power, speed and load angle). Oscillations are also propagated by AVR to excitation system values. Large system disturbances can cause significant changes to local parameters relevant for generator stability. Due to the linear nature of the voltage regulator system, a change in grid voltage or configuration leads to generator working in non optimal state. Even if the system returns to the same operating conditions, passing through nonlinear states during disturbance causes non optimal control, diminishing the stabilizing effects of the voltage regulator on the system. This can in turn cause another disturbance or, in case of larger units, disruption of the whole power grid. Limitations of the use of linear control systems call for an improved system which can provide optimal control of generator excitation system [13]- [15]. There is a lot of research effort directed toward design of excitation control system based on neural networks. However, very few of the demonstrated solutions have been implemented on real plants. The goal of the chapter is to develop a practical-to-use system. This is mainly achieved by separating the design on preparation and implementation phase. Almost all algorithm steps are executed in off-line mode, up to the final step of digital voltage regulator implementation. To implement the voltage regulator described in the chapter following equipment is needed: • data acquisition hardware and software • personal computer with Matlab and Simulink environment • automatic voltage regulator with neural network support (no training support is needed) Described controller is not adaptive. Adaptive critic design (ACD) algorithm is used for optimal control. Classic voltage regulators (AVR) equipped with PSS can be very effective in voltage control and power system stabilization if it is used in well configured transmission grid. In such conditions the generator stays in the linear operating range and the controller maintains optimal configuration on change of active and reactive power. Voltage regulator accuracy, balance of reactive power, stabilizing effect during large disturbances and active stabilization in normal operating conditions are achieved.

119

120

MATLAB Applications for the Practical Engineer

If the generator is connected to relatively weak transmission grid, with possible configuration changes during operation, use of neural network based control of excitation voltage leads to better results. 5.1. HDP algorithm implementation In order to achieve optimal control of synchronous generator connected to the power grid, it is necessary to select a minimum of two system variables – generator voltage and either output power or angular velocity. Selected variables must reflect system oscillations during distur‐ bances. It is also necessary to choose cost function that is minimized during training. System optimization is achieved by minimization of squared generator voltage governing error ΔU (t) = U G − U G and minimization of output power change ΔP obtained from: REF

U (t ) = ( K u1 × DU G (t ) + K u 2 × DU G (t  1) + K u 3 × DU G (t  2) ) + 2

( K P1 × DP(t ) + K P1 × DP(t  1) + K P1 × DP(t  2) )2

(27)

It is important to notice that the two demands are in conflict. Increase in gain in main control loop reduces static control error, but also reduces system damping and increases power output deviation ΔP . Use of cost function (27) ensures training process that increases system damping without sacrificing controller responsiveness. System identification is carried by training model neural network. Figure 8 shows training process of model neural network. Identification of nonlinear system is made by series-parallel structure. For ACD algorithm optimal control two output state space variables are chosen – generator voltage and output power, y = ΔU G ; ΔP . Vector y is used to form the cost function. Input vector is u (t ) = [ yr (t ) yr (t  1) yr (t  2) y (t  1) y (t  2) y (t  3)]

T

(28)

where yr = U GREF represents generator voltage setpoint,

uf – exciter control signal, y = ΔU G ΔP

T

is system output.

Model NN has 9 inputs and 2 outputs. Figure 9 shows simulation block diagram of identification procedure, implemented in Matlab and Simulink environment.

Dual Heuristic Neural Programming Controller for Synchronous Generator http://dx.doi.org/10.5772/58377

u (t)

[Ug] Mo d e l o f Syn ch ro n o u s g e n e ra to r

[Pg]

[Ug]

tm

[Ug]

1/z

[Pg]

1/z

TD

[Pg]

R N N Mo d e l

[Ug]

[Pg] V

[Ug]

Eva u la tio n

[Pg]

Figure 9. Model identification block diagram

A third order nonlinear mathematical model was used as a synchronous generator [20]. Synchronous generator is modeled as a single machine connected to an infinite bus. Although the model ignores initial magnetic flows, it proves to be a good approximation of an actual generator. A predictor-corrector method is used to solve differential equations. Selected mathematical model, used with well conditioned state space variables and described methods, achieves good numerical stability during simulation, allowing for relatively large sample time (TD=0.02s). Described model of synchronous generator is implemented in Clanguage S-function and used for further simulations. For model identification, recurrent neural network is used, with five neurons in hidden layer and two neurons in output layer. Hidden layer is implemented with full recursion and tansig activation function. Output layer uses linear activation function. Initial learning rate value η (21) is set to η = 105. It was reduced during training to η = 102. Training was completed in ~10000 steps. Results are shown in Figure 10. Signals U G NN(0) and U G NN(1) represent generator voltage at the start and at the end of RNN training. Figure 10a compares responses of control error of generator model and neural network. Voltage setpoint is generated as a pulse signal with a magnitude of ± 10%. In figure 10b, response to the change of generator power is shown. Figure 10c shows squared sum of RNN error during training. In the same experiment efficiency of the algorithm was verified. Training was stopped after 100000 steps. For the rest of the simulation, network output is calculated without the change of weight coefficient, thus eliminating the recency effect.

121

MATLAB Applications for the Practical Engineer

(a) 0.04

∆UG [p.u.]

∆ UGRE:F

∆UG NN(0)

∆ UG

∆ UG NN(1)

0.02

0 -0.02

0

1

2

3

4

5 t, s

6

7

8

9

10

(b) 0.06 P NN(0) P P NN(1)

∆P [p.u.]

0.04 0.02 0 -0.02 -0.04 0

1

2

3

4

5 t, s

6

7

8

9

10

(c) 0.01 0.008

Eror

122

0.006 0.004 0.002 0

1

2

3

4

5

6

7

8

9

10

11 4

k

x 10

Figure 10. Training results of model neural network

Figure 10.Training results of model neural network

Oncethe the model model neural neural network network is is trained, Once trained, itit is is possible possibletotocreate createsimulation simulationfor fortraining trainingcritic critic and in Figure Figure 11. 11. andaction actionneural neural networks. networks. Simulink Simulink model model of of the the simulation simulation is is shown shown in

[Pg]

[Ug]

x( t)

UG REF

TD SG

x(t+ 1)

TD

Ug

Pref

NN-Critic HDP

TD

Pg

e

Learn

tm th

Ru

RNiz

[Ug]

[dUdu]

SG

Ug

Referent

[Pg] User criteria U(t)

Pg

NN Model P

e

tm

[ dY1du]

[ dY2du]

Ug

[ dY1du] [ dY2du]

th

TD

NN action

[ dY1du]

[Ug]

Learn

Action learn

[ dY2du]

[dUdu]

[ Pg]

Figure 11. Simulink model for training action and critic networks

Figure 11.Simulink model for training action and critic networks

Dual Heuristic Neural Programming Controller for Synchronous Generator http://dx.doi.org/10.5772/58377

Simulation uses expressions (21), (22) and (23). Neural network training is implemented in Matlab S-functions, written in C programming language. Input of the S function is formed of system signals and desired value. S-function returns RNN ∂ yi (t) output signals and partial derivatives of outputs . Critic and action networks are ∂ uj (t)

)

(

trained in parallel.

Training results are shown in Figure 12. Figure 12a compares generator voltage response using conventional automatic voltage regulator (AVR) to neural network based controller trained using HDP algorithm. It can be seen that voltage rate of change was not degraded as a result of system damping. It is important to emphasize that generator voltage rate of change is key to system stability during large disturbances. Figure 12b compares active power output oscillation of conventional and HDP controller. π Figure 12c shows load angle, which is indicative of system stability δ < . System damping 2 is significantly increased by use of neural network controller.

(

)

D UG [p.u.]

(a) 0.04

RNN AVR

0.02

D UGREF

0 -0.02 -0.04 -0.06 0

1

2

3

4 (b)

0.15

5 t, s

6

7

8

9

RNN AVR

0.1

D P [p.u.]

10

0.05 0 -0.05 -0.1 0

1

2

3

4 (c)

5 t, s

6

7

8

9

10

d [rad]

0.7

0.6

0.5 AVR RNN

0.4 0

1

2

3

4

5 t, s

6

7

8

9

10

Figure 12. Response comparison of conventional AVR and HDP controller. a) generator voltage; b) generator active power; c) load angle

123

MATLAB Applications for the Practical Engineer

HDP algorithm is also verified on voltage control of synchronous generator from Simulink SimPower System block library (Synchronous Machine). Model network is designed with one hidden recursion layer with 5 neurons. Critic and action networks use the same structure, but with four neurons in hidden layer. HDP algorithm shown in Figure 11 is applied. Training of neural networks is performed in two steps– in first step action NN is trained to mimic classic voltage regulator, and in second step DHP algorithm is used to train both action and critic NN. It is important to get initial values of action NN in first step, as DHP needs network partial derivatives in training critic NN. Figure 13 shows training results. Figure 13a shows comparison of classic and NN controller responses. Generator voltage rate of change is almost the same for both controllers.

(b)

(a) 0.1

0.1

∆ UGREF ∆ UG AVR 0

0.05

∆ UG [p.u.]

∆ VG [p.u.]

∆ UGREF

∆ UG NN

0.05

∆ UG NN ∆ UG AVR

0

-0.05

-0.05

-0.1

-0.1 0

1

2

3

4

5 t, s (c)

6

7

8

9

0

10

0.15

1

2

3

4

5 t, s (d)

6

7

8

9

10

0.2 0.1

0.05

∆ P [p.u.]

∆ P [p.u.]

0.1

0 -0.05 -0.1

∆ UGREF

NN

0 -0.1

∆ UGREF

-0.2

AVR

NN

P AVR

-0.15 0

1

2

3

4

5 (e) t, s

6

7

8

9

0

10

1

2

3

4 (f)

0.8

5 t, s

6

7

8

9

10

1.2 1

∆ δ [rad]

0.6

∆ δ [rad]

124

0.4

∆ UGREF

0.2

NN AVR

0

0.8 0.6 0.4

∆ UGREF

0.2

NN AVR

0 0

1

2

3

4

5 t, s

6

7

8

9

10

0

1

2

3

4

5 t, s

6

7

8

9

10

Figure 13. Response comparison of conventional AVR and HDP controller a), b) generator voltage; c), d) generator active power; e), f) load angle

Figure 13c and 13e show comparison of output power and load angle change. It is obvious that the HDP neural network controller shows significant improvement over classic regulator.

Dual Heuristic Neural Programming Controller for Synchronous Generator http://dx.doi.org/10.5772/58377

Figure 13b shows comparison of voltage transient in the condition generator voltage reduced to 90% of nominal value. Rate of change of generator voltage is maintained even during such extreme operating conditions. Figure 13d shows conventional automatic voltage regulator (AVR) operating on the edge of stability, while the DHP controller is completely stable. Comparison of load angle values demonstrates significantly better behavior of DHP controller. On large load angles DHP controller maintains stability while conventional AVR provides low system damping, leading to instability and possible outages. Neural network HDP voltage controllers were developed using DHP algorithm. Experimental results showed that DHP algorithm is more efficient than ADC algorithm. Implementation of GDHP does not bring significant improvements over DHP. Application of developed S-functions enables simple use of described algorithms, as well as experimenting with different function parameters.

6. Real time windows target implementation DHP optimal control algorithm has also been implemented using Matlab Real Time Windows Target Toolbox (RTWT), based on [18] and [19]. Simulation hardware consists of desktop PC equipped with National Instruments data acquisition card NI DAQ 6024. Use of RTWT platform enables direct implementation of developed controllers on the real system. In Figure 14, a system consisting of synchronous generator and RTWT controller is shown. Signal acquisition is achieved using DAQ 6062E Input, with AD conversion of the input values. Line voltage and phase current signals are processed in “Clark” S-function, using Clarke transformation. Transformed values of generator voltage U G and generator power P are used as control feedbacks. Digital to analog conversion is performed by output module DAQ 6062 Output. All simulation blocks are part of a standard Simulink library except Clark transfor‐ mation and RNN S-functions.

ACD algorithm is implemented in two steps. Classic voltage controller is used in the first step (Figure 14). In the second step it is replaced with NN controller. During testing, generator active power and voltage were changed for a wide range of values. A pulse signal was superimposed to the generator voltage set-point value, resulting in large disturbances on the generator. To obtain measurement records RTWT Signal Visualization was used, with measured data saved in Workspace Simulink block yP. Recorded data was used to train neural network in offline mode. A RNN was used with one hidden layer containing 5 tansig neurons. Model NN had 9 inputs and 2 outputs. Results of NN training are displayed in Figure 15.

125

Figure 14.Automatic voltage control of synchronous generator using RTWT platform Figure 14. Automatic voltage control of synchronous generator using RTWT platform

D UG [p.u.]

During testing, generator active power and voltage were changed for a wide range of values. A pulse signal was superimposed to the generator voltage set-point value, resulting in large disturbances on the generator. To obtain measurement records RTWT Signal Visualization was used, with measured data saved in Workspace Simulink block yP. Recorded data was (a) 0.1 used to train neural network in offline mode. 0.05 A RNN was used with one hidden layer containing 5 tansig neurons. Model NN had 9 inputs 0 and 2 outputs. Resultsof NN training are displayed in Figure 15. -0.05 -0.1 -0.15

D UGREF

-0.2

D UG D UG NN

-0.25 -0.3

0

2

4

6

8

10 t, s

12

14

16

18

20

(b) 0.15 0.1 0.05

D PD UG [p.u.]

126

Clarke transformation. Transformed values of generator voltage U G and generator power P are used as control feedbacks. Digital to analog conversion is performed by output module DAQ 6062 Output. All simulation blocks are part of a standard Simulink library except Clark transformation and RNN S-functions. MATLAB Applications Practical Engineer ACD algorithmforisthe implemented in two steps. Classic voltage controller is used in the first step (Figure 14). In the second step it is replaced with NN controller.

0 -0.05 -0.1

D UGREF DP D P NN

-0.15 -0.2 0

2

4

6

8

10 t, s

12

14

16

Figure 15. Results of Model NN training a) generator voltage; b) generator active power

18

20

Dual Heuristic Neural Programming Controller for Synchronous Generator http://dx.doi.org/10.5772/58377

Obtained model NN was used to train critic and action RNN using procedure described in 3. The trained action NN was transferred to Simulink block (Figure 14) and prepared for use in real time. Using RTWT module system was run in real-time mode, with trained action RNN as a voltage controller.

(a) 0.9

U

GREF

UG [p.u.]

0.85

U AVR G

0.8

U NN G

0.75 0.7 0.65 0.6 0

1

2

3

4

5 t, s

6

7

8

9

(b)

55

P AVR P NN

50 45 P [kW]

10

40 35 30 25 20

0

1

2

3

4

5 t, s

6

7

8

9

10

(c) 100

δ AVR δ NN

δ[ °]

80 60 40 20 0

1

2

3

4

5 t, s

6

7

8

9

10

Figure 16. Classic controller and HDP NN controller – system responses: a) generator voltage; b) generator active pow‐ er; c) load angle

Simulation results are shown in Figure 16, comparing classic and DHP voltage controllers. Figure 16a represents comparison between generator voltage responses and shows satisfactory rate of generator voltage change and high control accuracy. In Figure 16b, generator active power outputs are compared. Considerable damping was achieved using DHP controller. Improved system stability during large disturbances is visible in Figure 16c, where the maximum load angle deviation is significantly lower for DHP controller. RTWT controller was implemented in sample-time mode with sample time set to T s = 2.5 ⋅ 10−4s. Control loop sample time is T D = 0.02s. Matlab Simulink RTWT Toolbox simplifies the transfer of developed procedures to the hardware platform.

127

128

MATLAB Applications for the Practical Engineer

7. Conclusion In the chapter, a complete process of development and implementation of neural network based controller is described. Optimal control is achieved by training neural network using data obtained from real system. Developed algorithms were tested on voltage control of synchronous generator connected to the power grid. Optimal control algorithm is implemented in Matlab and Simulink environment. Neural network training is implemented in C programming language Matlab S-function. S-function execution time was sufficiently low to achieve real time performance. Developed algorithm is tested using Matlab Real Time Windows Target Toolbox, using the same models with minimal modifications. Experiments performed on laboratory system show possibility of building optimal controller based on special network training procedure with no need for adaptive features. Obtained simulation results show advantages of neural networks based controller over classic controllers. Developed procedure can relatively easily be adapted for use on real systems.

Appendices Appendix A Matlab function for calculating hidden layer derivatives of outputs over weight coefficients of hidden layer. function[dadWR,a]=dadWR(WW1,WW2,p,a_1,dadWR_1,NW1,NW2); %function determines the partial derivative output recursive hidden%layers to weight coefficients matrix WW1% WW1 matrix of hidden layer% WW2 matrix of output layer% NW1 number of neurons in the hidden layer% NW2 length of the output vector neural networks% Ni length of the input vector neural networks% J target output value of the neural network% a output value of hidden layer neurons % p input vector of the hidden layer of recurrent neural network Ni=length(p)-NW1; a_temp=WW1*[p; 1]; a=tansig(a_temp); Ia=ones(NW1,1)-(a.*a); %value of derivative tansig activation function of neurons %y=WW2*a; %J=J-y; for l=1:NW1 for i=1:NW1 for j=1:Ni+NW1 dadar_temp=0; for k=1:NW1; dadar_temp=dadar_temp+WW1(l,Ni+k)*dadWR_1(k,(i-1)*Ni+j);

Dual Heuristic Neural Programming Controller for Synchronous Generator http://dx.doi.org/10.5772/58377

end dadWR(l,(i-1)*(Ni+NW1)+j)=Ia(i)*((l==i)*p(j)+dadar_temp); end, end, end, %end

function dadWR

Appendix B S-function for RNN training with one hidden layer. The first part defines sample time, number of inputs, outputs and neurons. Input is formed of values used for training. Output 1 is RNN output. Output 2 is partial derivative of output over input. Training is performed using Kalman filter. /* File : model.c * Abstract: * For more details about S-functions, see simulink/src/sfuntmpl_doc.c. * Copyright 1990-2002 The MathWorks, Inc. * $Revision: 1.12 $ */ /* ================model===================== * This sfunction is used for learning recurrent neural network (RNN) in programming environment Matlab Simulink * Learning RNN is implemented using the Kalman filter * This program supports RNN with one hidden layer * In the first part defines the time discretization, number of inputs, number of outputs and * the number of neurons in a single hidden layer * Simultaneously with learning RNN this function account the partial derivative of output per inputs * Port number of the s-function is equal to the number of inputs in RNN + number of outputs * Mato Miskovic Imotica */#define S_FUNCTION_NAME model #define S_FUNCTION_LEVEL 2 #include "simstruc.h" #include #include #define #define #define #define #define #define #define #define #define #define #define #define #define #define #define #define #define #define

TS 1 //sample time N_ULD 3 //number of inputs N_TD 1 //number of steps backwards (t-1), (t-3),(t-3), (t-N_TD) NUL (N_ULD*N_TD) NIZ 1 //number outputs W1 5 // number of neurons in the hidden layer W1_j (NUL+W1+1) // W2_j (W1+1) W1_ij (NUL+W1+1)*W1 W2_ij NIZ*W2_j N_W1 0 N_W2 N_W1+W1_ij SIZE_P (W1_ij+W2_ij) N_P N_W2+W2_ij NH_1 N_P+ SIZE_P*SIZE_P N_dadWR NH_1+W1 N_OUT N_dadWR+W1*W1_ij N_dnet N_OUT+NIZ

129

130

MATLAB Applications for the Practical Engineer

#define #define #define #define #define #define

N_ON N_dnet+NIZ*N_ULD NX N_ON+10 // ETTAMI 10.0 //final value of the learning rate ETTAQ .051 // Q ITER_UCI 20000 //number of steps to stop learning RNN NDISCSTATES NX

static void ttsig(real_T x, real_T *y)

//tansig

{ *y=2.0/(1.0+exp(-2.0*x))-1.0; } //AxB=C multiplying matrices in line form static void aputab_c(real_T *aa, real_T *bb, real_T *cc, int_T a_r, int_T a_c, int_T b_c){ int_T i, j, k; real_T priv; for(i=0; i