Welcome to the World of Modelling and Simulation

What is Modelling?

This blog is all about system dynamics modelling and simulation applied in the engineering field, especially mechanical, electrical, and ...

Application of the Taylor Series to Determine the Unknown Coefficients of a Polynomial

Let's say, we have a polynomial equation of the following form. 

y = ax^b + 5

Where, a, and b are the unknown coefficients of the polynomials. In this tutorial, we will discuss a method to determine the unknown coefficients using the Taylor series expansion assuming the data points are provided.

From the above polynomial equation, y is the dependent variable, and x is the independent variable. Let's assume that we know ten data points for the polynomial so that the following equations may be written based on the data points: 

y₁ = ax₁^b + 5

y2 = ax2^b + 5

y3 = ax3^b + 5

.
.
.

y10 = ax10^b + 5

This model seems to have a nonlinear dependence on the parameters a and b. Nonlinear regression can be used here to estimate these parameters. Nonlinear regression is based on determining the values of the parameters that minimize the sum of the squares of the residuals and the solution is preceded in an iterative manner. The Gauss-Newton method is one of the algorithms which minimize the sum of the squares of the residuals between nonlinear equation and data. A Taylor series expansion is used to express the nonlinear equation is an approximate linear form. After that, least square method is applied to estimate the parameters, which move towards the minimization of the residuals.

yi = f(xi; a,b,c) + ei;          i = 1, 2, 3, ... ... ... 10

Where, yi is a measured value of the dependent variable,  f(xi; a,b,c) is the nonlinear function of the independent variable, xi, with two unknown parameters a, b, and one given parameter c, which is 5 according to the given correlation. The last term ei is the random error. The model may be represented in a short form without the parameters for convenience as,

yi = f(xi) + ei                               (1)

For the two parameters case, a Taylor series expansion is written for the first two terms as,

Taylor series expansion for first two terms





Here, j is the initial guess, j+1 is the prediction, ∆a = aj+1 - aj and ∆b = bj+1 - bjEquation (2) is substituted to equation (1) which yields,

Taylor series expansion for first two terms




In matrix form, equation (3) may be written as,

Taylor series in matrix form




Where, Zj is the matrix of the partial derivatives of the function at the initial guess j.

Matrix of the partial derivatives


















The vector {D} contains the difference between the measurements and function values at 10 points.

Contains the difference between the measurements and function values











And, the vector {A} contains the changes in the parameters’ values as,

Contains the changes in the parameters






Applying least square method to equation (4) results the following normal equation,

Applying the least square method




Equation (5) is implemented to solve for {∆A} to estimate the improved values of the parameters a and b,

aj+1 = aj + ∆a

bj+1 = bj + ∆b

Thus, this procedure is continued until the solution converges.



#TaylorSeries #Polynomial #NonlinearRegression #LeastSquare #Matlab #Blog #Blogger

Discussion of the Secant Method to Solve an Equation

Let's consider the following equation:

eˣ = 2x² + 1

We may write the above expression as,

f(x) = eˣ - 2x² - 1

To find the value of x or roots of the equation, we may apply the Secant method. In this method, the derivative may be approximated by a backward finite divided difference approach. So, the function may be expressed as,
Secant method expression
This method is convenient when evaluating derivative for some function is difficult in Newton-Raphson method. Now, the above equation may be substituted in Newton-Raphson’s formula to get the following algorithm.
Secant formula
In contrast to Newton-Raphson method, two initial assumptions are made for the two unknowns’ xi-1 and xi.
For the proper initial guess for x, graphical method is applied to visualize the characteristics of the function. The MATLAB commands generate the plot of the function f(x). This function is plotted with respect to the values of x from -10 to +10. The following figure shows the function plot. From the plot, the function changes sign when x is between 2 to 3. So, the root lies in this interval. And, in the program developed for Secant method, two initial guesses for xi-1 and xi are 2 and 3 respectively.

Plot of the function f(x) to identify the interval where root of the function lies














However, more close observation reveals two roots, which are shown in the following figure. In this time, the range of x is decreased significantly for precise investigation between -2 to +2 as this interval seems unclear in the previous figure for the presence of any roots. 

Shows two more roots which lie in between 0.1 to 1, and at origin















So, the above figure shows the existence of other two roots, which lie in between 0.1 to 1 and at the origin. The origin is obvious for this case, and there is no need for any further estimation by numerical methods. But, for the interval 0.1 to 1, the accurate and precise approximation is necessary. For this, the two initial guesses are set to 0.1 and 1 for xi-1 and xi accordingly.
The following Secant formula is implemented to approximate the two roots lie in the intervals 2 to 3 and 0.1 to 1. The results are shown after the program.


Secant Formula Implementation by a MATLAB Program

% Definition of a function "secant" to solve equation by Secant Method
 
function [x,ea] = secant(X,X0,etol)
 
format long;
 
% Input: X,X0, etol=Tolerance definition for error
% Output: x, ea=Calculated error in loop
% Iterative Calculations
 
while (1)
    
    % Implementation of Secant Algorithm 
          
solution = X-((exp(X)-(2*X^2)-1)*(X0-X))/((exp(X0)-(2*X0^2)-1)-(exp(X)-(2*X^2)-1));     
    solutionprevious=X;
    X0=X;
    X=solution;
       
           
if (solution-solutionprevious)~=0        % Approximate percent relative error
                                 
   ea=abs((solution - solutionprevious)/solution)*100;
       
end
 
if ea<=etol,                  
   
    break,                
end
 
x = solution;
 
disp(x);
 
end
 
% Display of output parameters
 
disp(x);                      
 
disp(ea);
 
%disp(solution-(0.1*abs(ea)*solution));
end


Outputs:

>> secant (3, 2, 1e-5)
% Iterations

   2.597424571529533
   2.796706966294471
   2.858723109017894
   2.841817367110900
   2.842657964680533
   2.842673428005184
   2.842673428005184

>> secant(1, 0.1, 1e-5)
% Iterations

   0.308929151717719
   0.570046187158243
   1.157486960694727
   0.682985699401547
   0.723807551583747
   0.742065590882227
   0.740827136083267
   0.740850398734311
   0.740850398734311


So, the solutions for the equation are 2.842673428005184, 0.740850398734311 and 0.



#SecantMethod #RootsofEquations #Matlab #NumericalMethod #Blog #Blogger

Modelling and Simulation of an Electric Vehicle in SIMULINK

In this tutorial, I am going to present you a very simple way to design and control your very first electric vehicle in SIMULINK. You need to have MATLAB version 2019/2020 with full features. Most of the features that I will be using are from Simulink and Simscape libraries. Simscape is a very powerful tool that Matlab has integrated lately. This tutorial will be straightforward - showing you the diagrams and explaining afterwards. At first, let's begin with a simple question, how an overall electric vehicle system looks like? Below is the image of the electric vehicle system. 

Complete Simulink-Simscape block diagram of an electric vehicle system












As it appears, this is a Simulink block diagram, which contains many blocks from the library browser. Let's see what are the fundamental components of an electric vehicle system. From the diagram, there are five main blocks:

1. Vehicle Body 
2. DC Motor
3. DC Motor Driver
4. Battery
5. Controller

1. Vehicle Body

In this modelling, we can consider some basic properties of a vehicle. Below, we see the model of the vehicle body, which has four tires, a simple gear, and vehicle property block. When you click on each block, you will find the respective properties that you need to choose or enter independently. Just remember, these blocks represent the governing dynamics that Matlab has already prepared for us in the Simulink blocks. You do not need to spend hours for deriving system equations. Rather, here, we ought to be paying attention while connecting the blocks and providing appropriate parameters.

Representing vehicle body by Simscape blocks
















When you click on the vehicle body block, it will pop-up a new window as below. For this tutorial, I just chose the default parameters, but you can certainly change these for your specific problem.

Simscape vehicle body block parameters
















The following is the tire property block.

Simscape vehicle tire model block


















This is how a simple gear block appears when you double click the block.

Simscape simple gear box model block
















2. DC Motor

The purpose of the DC motor in electric vehicle is to drive the tires/wheels. It should be noted that from the Simscape library, you can choose different types of motors. The one side of the DC motor block is connected to the vehicle body blocks, and the other side is connected to the DC motor driver, which controls and rotates the motor. The following image shows the DC motor block in Simulink.

DC motor block properties


3. DC Motor Driver

To precisely control a DC motor, we are required to design a driver for that. There are numerous ways you may develop the motor controller. I used an H-Bridge block from Simscape library along with the PWM voltage controller block. In the following image, you see that there are some other blocks, such as, current and voltage sensors, electrical reference etc. 

DC motor driver system in Simulink - Simscape















4. Battery

For any electric vehicle, battery is a very important part. Still, in industries, people are struggling to find a better solution for the battery that can run an electric vehicle with more mileage before it drains fully. Here, I have chosen a simple battery model from Simscape library where you can obviously tune the parameters for your design requirement. To implement the battery and connect it to the rest of the system, a controlled current source block is added. As before, we need to have an electrical reference for the current source block, which is depicted in the subsequent figure.

Simscape battery system blocks

5. Controller

In newer versions Matlab with full features, you will find a block named 'Longitudinal Driver' that may act as a controller for the electric vehicle system. Below, we see the block parameters. I am using all the default properties for this block. You can decide the type of the controllers at first, then shift type, units, output gear signal, and so forth. There is a whole lot option to tune this block parameter for your specific goals.

Showing longitudinal driver block properties





















We begin this tutorial with the final version of the Simulink system model of an electric vehicle. There were some blocks that were masked by some images, below you can see the original version of the complete system diagram for your first virtual electric vehicle model. As you can see, in this model, most of the blocks are from the Simscape library, and you cannot connect them directly to any other Simulink blocks. But, with a converter, you can easily make the connections. In Simscape, there are two blocks - PS-Simulink and Simulink-PS converter - which are used to convert physical signals from Simscape block into Simulink compatible signal and vice versa.

Complete Simulink-Simscape block diagram of an electric vehicle system













Some Results

Everything is done. Now, if you would like to see the responses of different parts of the electric vehicle, you can do so. You need to add the oscilloscope blocks from the Simulink library to those parts. 

Here is the source signal that is fed into the system to run the vehicle.

Signal generated from the source
















This is the signal that is an output from the longitudinal driver block.

Signal from the longitudinal driver block






















This oscilloscope shows the current from a current source block.

Shows the current from a current source block.



















That's all. We have designed a virtual model for an electric car. Certainly, this is not a Tesla as the picture of the vehicle body block shows in the first image. Nevertheless, we can at least have a little sense of accomplishment that a simple model may even be designed in our laptops. 



#ElectricVehicle #Simulink #Simscape #Modelling #Simulation #DCmotor #Blog #Blogger

Medical Image Segmentation with 3D Slicer: A Beginner's Guide

3D Slicer is an open source software that is widely used for image processing, visualization, filtering, and so forth. In this article, I would like to discuss its effectiveness and application in medical image segmentation. This software is completely FREE, and there is a continuous development of many add-ins/extensions for this platform. You may download the software from 3D Slicer Webpage. After downloading, the installation procedure is straightforward, you may also follow the instructions mentioned in their webpage.

Medical image segmentation to extract the size or volume of an organ or complex airways/channels from computed tomography (CT) or micro-computed tomography (𝛍CT) is very interesting and has been playing a crucial part in biomedical engineering. For example, human nasal cavities or airways have such a complex formation that from the CT scans, we are unable to extract the volume. However, with the power of modern image processing tools, we can manually trace the nasal airways path slice by slice and finally obtain the volume of the cavities.

In this tutorial, I am going to highlight the 3D Slicer's interface, and how we can do the manual/automatic segmentation for the purpose of extracting the volume or shape information. When you launch the software, you will see the following welcome window. The window has four sectional views. The three black windows show the transverse (or, axial), sagittal (or, longitudinal), and coronal (or, frontal) views of the CT images. The other blue colored section highlights the three dimensional reconstruction of the segmented images.

3D Slicer Interface after Launching
















Medical images are often processed as .DICOM file, which can be easily opened in 3D Slicer. In the welcome window, there are options to load different format of data, but Digital Imaging and Communications in Medicine (DICOM) images are smoothly handled by the 3D Slicer. In the following, I have opened a CT scan of a human brain with a tumor, which is available in 3D Slicer's sample data for demonstration. Our task here is to find the volume of the tumor quantitatively from the CT scan images so that physicians can have more in-depth information regarding the size and growth of the tumor. They can make their next decisions confidently based on these data.

Loading the CT images in 3D Slicer
















To do so, at first, we need to open the segment editor, which is located at the top of the window where you will see an icon says 'All Module'. Click on it and then, select the segment editor as highlighted below.

Segment Editor of the 3D Slicer
















Once you select the segment editor, a new small window will pop-up just like the screenshot shown beneath. You will find some features in that new little window, such as, 'ADD', 'Show 3D' etc.

Adding a New Segment to Start the Segmentation
















You need to choose the option 'Add (+)' and then, a new segment will be created named as 'Segment_1' (see below). From the screenshot, you see that a new segment is ready to carry on for further operation.

A New Segment Created in the Property Editor
















Now, if you double click on the 'Segment_1', you will see the following options will appear to select the segmentation nature. For example, if we segment a bone, then choose 'bone' in that window; or if the segmentation is related to the brain or tooth, select the respective 'brain' and 'tooth' segmentation features. In our example, it is a tumor that we are going to segment. A tumor is a mass of abnormally growth of tissues, so we can choose the segmentation media as 'tissue'.

Selection of the Type or Nature of Segmentation
















So, we have decided how we are going to segment the tumor, now it's time to move forward using the tools to do the segmentation. There are two types of segmentations: automatic and manual. To understand the difference between the automatic and manual segmentation, let me first explain how the CT scans are formatted or processed for image analysis. There are many single images that form the CT scan of an organ, which alternatively implies that a CT scan itself an image that has subsets of many images. So, when we consider automatic segmentation, it is a process of tracing or narrowing a particular section from the CT scan that is done by the software, and manual segmentation is that you need to paint or trace that particular region by hand for all the individual images that eventually form the CT scan. As you see below, the red circled area is the tumor that we need to separate out and get the volume of it. We also see that there are some options, like 'Paint', 'Draw', 'Level Tracing' to trace the tumor. Here, I used the 'Level Tracing' option for that, traced few images or slices and then used the 'Fill between slices' option as an automatic segmentation approach.

Use of the Segmentation Tools

Once the segmentation is done, then we can estimate the volume of the tumor. For that, select again 'All Modules', and from the dropdown menu, click on the 'Quantification', and then hit 'Segment statistics'.


How to Extract Volume Information from 3D Slicer
















After that, at the bottom of your window, you will find the volume of the segmented tumor. From the screenshot below, we see the volume of the tumor is approximately 346.99 mm³. Now, remember, this is not an accurate or quantitatively perfect representation of the tumor. However, it gives a reasonable view and idea about the shape and size of the tumor. Below, we can see the looks of the tumor in a three dimensional space.

The 3D Shape Extracted from Segmentation with Volume Statistics



#ComputedTomography #3DSlicer #ImageSegmentation #ImageProcessing #3DModelling #Blog #Blogger

A Taylor Series Expansion to Approximate a Function near a Point

The Taylor series expansion up to seventh terms for approximation of near point x =1 is,
Taylor series expansion up to seventh terms





Letting starting point xₒ is at 0 and the step size (xi - xₒ) is 1. Now, the value of eˣ at x = 1 is to be calculated by Taylor series expansion on the basis of the starting point and its derivatives at x = 0 with step size 1. A MATLAB program is written to evaluate this value from zero order approximation to sixth order approximation and show the improvement of accuracy and minimization of truncation error as well.
The true value of is 2.718281828459046. The truncation error is the difference between true and approximate values. The following MATLAB program describes each steps of the calculation. 

Taylor Series Implementation by a MATLAB Program

% Taylor Series expansion to approximate exp(x) near point x=1
 
function e = taylorseries(x0,x1,h)
 
e0=exp(x0);              % calculates exp(x) value when x=0
 
derivative1=1;           % all the six derivatives calculated at x=0
derivative2=1;
derivative3=1;
derivative4=1;
derivative5=1;
derivative6=1;
 
% zero order approximation
 
        e=e0;   
        error=exp(x1)-e;      
        disp (e);        % display of exp(x) and error at command prompt          
        disp(error);
        
% first order approximation        
            
        e=e0+(derivative1*h);       
        error=exp(x1)-e;       
        disp (e);       
        disp(error);
             
% second order approximation
 
        e=e0+(derivative1*h)+(derivative2*(h^2/factorial(2)));        
        error=exp(x1)-e;        
        disp (e);        
        disp(error);
        
% third order approximation
                     e=e0+(derivative1*h)+(derivative2*(h^2/factorial(2)))+(derivative3*(h^3/factorial(3)));        
        error=exp(x1)-e;        
        disp (e);        
        disp(error);
        
% fourth order approximation        
                  e=e0+(derivative1*h)+(derivative2*(h^2/factorial(2)))+(derivative3*(h^3/factorial(3)))+(derivative4*(h^4/factorial(4)));       
        error=exp(x1)-e;        
        disp (e);        
        disp(error);
        
% fifth order approximation                      
        e=e0+(derivative1*h)+(derivative2*(h^2/factorial(2)))+(derivative3*(h^3/factorial(3)))+(derivative4*(h^4/factorial(4)))+(derivative5*(h^5/factorial(5)));        
        error=exp(x1)-e;        
        disp (e);        
        disp(error);
        
% sixth order approximation
        e=e0+(derivative1*h)+(derivative2*(h^2/factorial(2)))+(derivative3*(h^3/factorial(3)))+(derivative4*(h^4/factorial(4)))+(derivative5*(h^5/factorial(5)))+(derivative6*(h^6/factorial(6)));        
        error=exp(x1)-e;        
        disp (e);        
        disp(error);
end


Table 1: Outputs from the Above Program to Approximate at x=1

Approximation of an exponential function by Taylor series










The above table tells us as the order of the Taylor series increases the truncation error also decreases. However, after sixth and higher order approximations, improvement of the accuracy is not much significant. So, we can stop the approximation to first seven terms of the Taylor series and the approximate value of ex near point x = 1 is 2.718.

Finding Roots of a Equation by Newton-Raphson Method

In this tutorial, we are interested to the roots of the following function.

f(x) = sin x + 2 - 0.1 x

To find the value of x or roots of the equation, Newton-Raphson method is applied. For the proper initial guess for x, graphical method is applied to visualize the characteristics of this nonlinear function. This graphical technique also helps to find multiple roots for the unknown function. The following MATLAB commands generate the plot of the function, f(x). This function is plotted with respect to the values of x from -50 to +50.

>> X=-100:100;
>> plot(X,(sin(X)+2-(0.1*X)))

From the following figure 1, it is clear that the equation has seven roots and plotting the function helps to get the roots. The function changes its sign from, x = 10 - 11, 11 -12, 16 -17, 18 - 19, 21 - 22, 25 - 26, and 27 - 28. And, the roots of the equation must lie between these intervals. Now, to get good approximations of the roots, a Newton-Raphson algorithm is written in MATLAB.

Plot of the function f(x) to identify the intervals where roots of the function lie















Figure 1: Plot of the function f(x) to identify the intervals where roots of the function lie.

The program is developed by MATLAB user defined function called “sinusoidal”. The initial assumptions for approximating seven roots of x are 10.5, 11.5, 16.5, 18.5, 21.5, 25.5 and 27.5 respectively. The codes are given below. It is straightforward to run the program seven times as we have seven initial guesses for the corresponding roots.

Newton-Raphson Algorithm by MATLAB

% Definition of a function "sinusoidal" to solve equation by Newton-Raphson
function [x,ea] = sinusoidal(X, etol)
format long;
 
% Input: X=21.5, etol=Tolerance definition for error
% Output: x=Solution of equation, ea=Calculated error in loop
 
% Program Initialization
 
solution = 21.5;
 
% Iterative Calculations
 
while (1)
    solutionprevious = solution;

% Implementation of Newton-Raphson Method 
           
    solution = X - (sin(X)+2-(0.1*X))/(cos(X)-0.1);               
    X=solution;       
           
if solution~=0                         % Approximate percent relative error
                                 
    ea=abs((solution - solutionprevious)/solution)*100;      
end
 
if ea<=etol,                  % Condition to meet specified error tolerance
   
    break,                
end
 
x = solution;
 
disp(x);
 
end

% Display of output parameters
 
disp(x);                      
disp(ea);
 
end

Program Outputs: 

>> sinusoidal (10.5,1e-10)
ans =
  10.636782424281707
>> sinusoidal(11.5, 1e-10)
ans =
  11.562055865585652
>> sinusoidal (16.5,1e-10)
ans =
  16.107752970689450
>> sinusoidal (18.5,1e-10)
ans =
  18.721338781142521
>> sinusoidal (21.5,1e-10)
ans =
  21.809224297055675
>> sinusoidal (25.5,1e-10)
ans =
  25.744697403517950
>> sinusoidal (27.5,1e-10)
ans =
  27.435908953388580

So, the seven roots from the program are 10.637, 11.562, 16.108, 18.721, 21.809, 25.745 and 27.436. In all the cases, the error tolerance limit is kept extremely small (0.0000000001) to get more accurate and precise results.


#Blog #Blogger #NewtonRaphson #Matlab #Sinusoidal