%***************************************************** % This file contains code for solving some % exercises, design problems, and examples in the book: % % Kevin M. Passino and Stephen Yurkovich, Fuzzy Control, % Addison Wesley Longman, Menlo Park, CA, 1998. % % The authors of each program are given at the % beginning of each program. % % Software Version: 7/6/98 % This is the last date that any program contained in % this file was changed; the latest version of this % file is kept at the web and ftp sites for the book: % http://www.awl.com/cseng/ % ftp://ftp.aw.com/cseng/authors/passino/fc % (use anonymous ftp). % %***************************************************** % Instructions: % % (1) Below, there are many different programs, where % each one is used for solving a particular problem or % example in the book. They were all grouped into this % one file for ease of transfer. Hence, to use the % computer programs below you will want to copy the % relvant program out of this file and give it an % appropriate name with an appropriate extension (e.g., % fuzzycontroller.c for a C program or fmrlc.m for % a Matlab program). % % (2) The programs are separated by markers: %***************************************************** >>--->>CUT-HERE>>--------------------------------- %***************************************************** % so that you can easily search for these and copy % the programs that you need. Note that each program % is labeled with the particular exercise or design % problem that it solves. Note that for brevity % EX.Y is used to denote Exercise X.Y and DPX.Y is % used to denote Design Problem X.Y. % % (3) All the programs are written in Matlab except % the ones written in C for DP2.1 and E6.1 (which is % programmed in both Matlab and C). % % (4) Be careful to study the relevant material in % the text and read the homework problem carefully % so that you know how to run/modify the program. % % (5) You are encouraged to modify the programs for % your own purposes. Note that even though many % of the programs are for relatively simple problems % the basic structure of the code will not change % too much when you extend/update the code to make % it suitable for solving complex real-world problems. % %***************************************************** % Instructional Material Disclaimer: The programs given % below are for their instructional value. They have % been tested with care but are not guaranteed for any % particular purpose. Neither the publisher or the % authors offer any warranties or representations, % nor do they accept any liabilities with respect to % the programs. %***************************************************** >>--->>CUT-HERE>>--------------------------------- %***************************************************** /* DP2.1(a) */ /* Program: fc.c Written by: Scott Brown, 5/1/95 2 input fuzzy controller to control inverted pendulum system. Controller has 5 membship functions for each input and 5 membership functions for the output. Center-of-gravity is used for defuzzification. */ #include #include #include #define MAX(A,B) ((A) > (B) ? (A) : (B)) #define MIN(A,B) ((A) < (B) ? (A) : (B)) #define PI 3.14159265359 /************************************************************************************/ typedef struct in_mem { double width; /* Input membership function width (1/2 of triangle base). */ double *center; /* Center of each input membership function. */ double *dom; /* Degree of membership for each membership function. */ } IN_MEM; typedef struct out_mem { double width; /* Output membership function width (1/2 of triangle base). */ double *center; /* Center of each output membership function. */ } OUT_MEM; typedef struct fuz_sys { IN_MEM *emem; /* Groups all fuzzy system parameters in a single variable. */ IN_MEM *edotmem; OUT_MEM *outmem; } FUZ_SYS; /************************************************************************************/ /* Function Prototypes: */ void fuzzy_init(FUZ_SYS *fuzzy_system); void fuzzy_free(FUZ_SYS *fuzzy_system); double fuzzy_control(double e, double edot, FUZ_SYS *fuzzy_system); void fuzzyify(double u, IN_MEM *mem); double leftall(double u, double w, double c); double rightall(double u, double w, double c); double triangle(double u, double w, double c); void match(const IN_MEM *emem, const IN_MEM *edotmem, int *pos); double inf_defuzz(IN_MEM *emem, IN_MEM *edotmem, OUT_MEM *outmem, int *pos); /************************************************************************************/ void fuzzy_init(FUZ_SYS *fuzzy_system) { /* Define the input and output membership functions. */ int i; /* Allocate memory for membership functions. */ if (!(fuzzy_system->emem = (IN_MEM *) malloc(sizeof(IN_MEM)))) { printf("Error allocating memory.\n"); exit(1); } if (!(fuzzy_system->edotmem = (IN_MEM *) malloc(sizeof(IN_MEM)))) { printf("Error allocating memory.\n"); exit(1); } if (!(fuzzy_system->outmem = (OUT_MEM *) malloc(sizeof(OUT_MEM)))) { printf("Error allocating memory.\n"); exit(1); } if (!(fuzzy_system->emem->center = (double *) malloc(5*sizeof(double)))) { printf("Error allocating memory.\n"); exit(1); } if (!(fuzzy_system->emem->dom = (double *) malloc(5*sizeof(double)))) { printf("Error allocating memory.\n"); exit(1); } if (!(fuzzy_system->edotmem->center = (double *) malloc(5*sizeof(double)))) { printf("Error allocating memory.\n"); exit(1); } if (!(fuzzy_system->edotmem->dom = (double *) malloc(5*sizeof(double)))) { printf("Error allocating memory.\n"); exit(1); } if (!(fuzzy_system->outmem->center = (double *) malloc(5*sizeof(double)))) { printf("Error allocating memory.\n"); exit(1); } /* Initialize for inverted pendulum. */ fuzzy_system->emem->width = PI/4.0; /* Width defined to be 1/2 of triangle base. */ fuzzy_system->edotmem->width = PI/8.0; fuzzy_system->outmem->width = 10; for (i=0; i<5; i++) { fuzzy_system->emem->center[i] = (-PI/2.0 + i*PI/4.0); fuzzy_system->edotmem->center[i] = (-PI/4.0 + i*PI/8.0); fuzzy_system->outmem->center[i] = (-20.0 + i*10.0); } } /************************************************************************************/ void fuzzy_free(FUZ_SYS *fuzzy_system) { /* Free memory allocated in fuzzy_init(). */ free(fuzzy_system->emem->center); free(fuzzy_system->emem->dom); free(fuzzy_system->edotmem->center); free(fuzzy_system->edotmem->dom); free(fuzzy_system->outmem->center); free(fuzzy_system->emem); free(fuzzy_system->edotmem); free(fuzzy_system->outmem); } /************************************************************************************/ double fuzzy_control(double e, double edot, FUZ_SYS *fuzzy_system) { /* Given crisp inputs e and edot, determine the crisp output u. */ int pos[2]; fuzzyify(e, fuzzy_system->emem); fuzzyify(edot, fuzzy_system->edotmem); match(fuzzy_system->emem, fuzzy_system->edotmem, pos); return inf_defuzz(fuzzy_system->emem, fuzzy_system->edotmem, fuzzy_system->outmem, pos); } /************************************************************************************/ void fuzzyify(double u, IN_MEM *mem) { /* Fuzzify the input u by determining the degree of membership for each membership function in mem. Assumes 5 membership functions, with first and last membership functions leftall and rightall respectively. Other membership functions are triangular. */ int i; mem->dom[0] = leftall(u, mem->width, mem->center[0]); for (i=1; i<4; i++) mem->dom[i] = triangle(u, mem->width, mem->center[i]); mem->dom[4] = rightall(u, mem->width, mem->center[4]); } /************************************************************************************/ double leftall(double u, double w, double c) /* Determine degree of membership for a leftall membership function. NOTE: u is input, c is mem. fun. center, and w is mem. fun. width. */ { if (u < c) return 1.0; else return MAX(0,(1-(u-c)/w)); } /************************************************************************************/ double rightall(double u, double w, double c) /* Determine degree of membership for a RIGHTALL membership function NOTE: u is input, c is mem. fun. center, and w is mem. fun. width. */ { if (u >= c) return 1.0; else return MAX(0,(1-(c-u)/w)); } /************************************************************************************/ double triangle(double u, double w, double c) /* Determine degree of membership for a TRIANGLE membership function NOTE: u is input, c is mem. fun. center, and w is mem. fun. width. */ { if (u >= c) return MAX(0,(1-(u-c)/w)); else return MAX(0,(1-(c-u)/w)); } /************************************************************************************/ void match(const IN_MEM *emem, const IN_MEM *edotmem, int *pos) { /* For each universe of discourse, determine the index of the first membership function with a non-zero degree (i.e. match the rules to the current inputs to find which rules are on). These indices are used to determine which four rules to evaluate. (NOTE: A 2 input sytem with no more than 50% overlap for input membership functions only requires the evaluation of at most 4 rules.) */ int i; for (i=0; i<5; i++) { if(emem->dom[i] != 0) { pos[0] = i; break; } } for (i=0; i<5; i++) { if(edotmem->dom[i] != 0) { pos[1] = i; break; } } } /************************************************************************************/ double inf_defuzz(IN_MEM *emem, IN_MEM *edotmem, OUT_MEM *outmem, int *pos) { /* We use the degrees of membership found in the function match() to form the implied fuzzy sets. The correct output membership function for each rule is determined by adding (and saturating) a shifted version of the input membership function indices (this implements the typical pattern of linguistic-numeric indices in the body of the table of rules). In this way we compute the rule-base at every step, rather than storing the rule-base in a table. Defuzzification is also performed using the center-of-gravity method. A crisp output is returned. */ double outdom, area, Atot = 0, WAtot = 0; int i, j, out_index; for(i=0; i<2; i++) { for(j=0; j<2; j++) { if ( ((pos[0]+i)<5) && ((pos[1]+j)<5)) { /* Check that bounds are not exceeded. */ outdom = 0; /* Shift indices left. */ out_index = ((pos[0]+i)-2) + ((pos[1]+j)-2); /* Saturate */ if (out_index < -2) out_index = -2; else if (out_index > 2) out_index = 2; /* Shift indices right.*/ out_index += 2; /* Determine the certainty of the premise */ outdom = MIN((emem->dom[pos[0]+i]), (edotmem->dom[pos[1]+j])); /* Defuzzify */ area = 2*outmem->width*(outdom - (outdom*outdom)/2); Atot += area; WAtot += area*outmem->center[out_index]; } } } /* Return the crisp value. Minus sign required to give correct output for pendulum system! Note that this minus sign actually ensures that the table of indices works out as shown in class. */ return -(WAtot/Atot); } /************************************************************************************/ #if 0 /* Set to zero to call fuzzy_control() from another program. */ /* Test I/O behavior of fuzzy controller. */ void main(void) { /* Test for input given below. Output is -6.818182. */ double e, edot, u; FUZ_SYS fuzzy_system; /* Crips inputs. */ e = 0; edot = PI/8.0 - PI/32; fuzzy_init(&fuzzy_system); u = fuzzy_control(e,edot,&fuzzy_system); fuzzy_free(&fuzzy_system); printf("e = %f, edot = %f, u = %f\n",e,edot,u); } #endif /************************************************************************************/ %***************************************************** >>--->>CUT-HERE>>--------------------------------- %***************************************************** /* DP2.1(a) */ /* Program: rk4.c Written by: Scott Brown, 5/1/95 Fourth Order Runge-Kutta Algorithm, for use with single-input systems. Function rk4(double *t, double *y, double u, double h, int n) is called with the 5 arguments indicated. Note that the first two arguments are pointers. Also note that no values are returned. Derivatives must first be defined in function deriv. Time is updated within rk4 call, so should not be updated by programmer in calling loop. The arguments for rk4() are defined below: double *t : pointer to time double *y : pointer to states double u : system input (for discrete simulations) double h : integration step size int n : system order The parameter u is useful for simulating discrete-time systems or physical systems in a laboratory. For continuous system simulation the input should be calculated within the deviv() function. */ /*************************************************************************** *********/ #include #include #include /*************************************************************************** *********/ /* Data Types */ typedef struct rk_struct { double *f1; /* Pointers to the Runge-Kutta coefficient arrays. */ double *f2; double *f3; double *f4; double *arg2; /* Pointer to a temporary vector used in several functions. */ } RK_STRUCT; /*************************************************************************** *********/ /* Application dependent area; change as necessary. */ #include "fuzzy.c" /* Include fuzzy functions. */ /* Global Variables */ FUZ_SYS fuzzy_system; /* Define fuzzy_system as a global variable. */ double input; /* System input; used to store input to file. */ /*************************************************************************** *********/ /* Function Prototypes */ void rk4(RK_STRUCT *rks, double *t, double *y, double u, double h, int n); void deriv(double t, const double *y, double u, double *ydot); void F1(double t, const double *y, double u, double h, int n, double *f1); void F2(double t, const double *y, double u, double h, int n, const double *f1, double *f2, double *arg2); void F3(double t, const double *y, double u, double h, int n, const double *f2, double *f3, double *arg2); void F4(double t, const double *y, double u, double h, int n, const double *f3, double *f4, double *arg2); void rk_init(RK_STRUCT *rks, int n); void rk_free(RK_STRUCT *rks); /*************************************************************************** *********/ /* Implementation */ void deriv(double t, const double *y, double u, double *ydot) { /* This is where the derivative is defined: xdot = f(t,x,u). */ /* Example: xdot = A*x + b*u (time invariant system) ydot[0] = -2*y[0] + -u; | -2 0 0 | | -1 | ydot[1] = -3*y[1] + 2*u; A=| 0 -3 0 |, B=| 2 | ydot[2] = -4*y[2] + u; | 0 0 -4 | | 1 | */ /******* INSERT YOUR DERIVATIVE DEFINITION HERE: ********/ double e,edot; /* Determine the system input u. Because the system is continuous, the input is calculated within the deriv() function. For discrete time systems, the input can be calculated and passed to deriv() via the u paramater. */ e = -y[0]; edot = -y[1]; u = fuzzy_control(e,edot,&fuzzy_system); input = u; /* Store u in input for plotting. */ /* Inverted pendulum dynamics. */ ydot[0] = y[1]; ydot[1] = (9.8*sin(y[0]) + cos(y[0])*(-y[2] - 0.25*pow(y[1],2.0)*sin(y[0]))/1.5)/(.5*(4.0/3.0 - 1.0/3.0*pow(cos(y[0]),2.0))); ydot[2] = 100.0*u -100.0*y[2]; /* Actuator dynamics */ } /*************************************************************************** *********/ void rk4(RK_STRUCT *rks, double *t, double *y, double u, double h, int n) { int i; /* Determine Runge-Kutta Coefficients. */ F1(*t,y,u,h,n,rks->f1); F2(*t,y,u,h,n,rks->f1,rks->f2,rks->arg2); F3(*t,y,u,h,n,rks->f2,rks->f3,rks->arg2); F4(*t,y,u,h,n,rks->f3,rks->f4,rks->arg2); /* Update time and output values. */ *t += h; for (i = 0; i < n; i++) { y[i] += (1.0/6.0)*(rks->f1[i] + 2*(rks->f2[i]) + 2*(rks->f3[i]) + rks->f4[i]); } } /*************************************************************************** *********/ void F1(double t, const double *y, double u, double h, int n, double *f1) { int i; deriv(t,y,u,f1); for (i = 0; i < n; i++) { f1[i] = h*f1[i]; } } /*************************************************************************** *********/ void F2(double t, const double *y, double u, double h, int n, const double *f1, double *f2, double *arg2) { int i; for (i = 0; i < n; i++) { arg2[i] = y[i] + 0.5*f1[i]; } deriv((t+0.5*h), arg2, u, f2); for (i = 0; i < n; i++) { f2[i] = h*f2[i]; } } /*************************************************************************** *********/ void F3(double t, const double *y, double u, double h, int n, const double *f2, double *f3, double *arg2) { int i; for (i = 0; i < n; i++) { arg2[i] = y[i] + 0.5*f2[i]; } deriv((t+0.5*h), arg2, u, f3); for (i = 0; i < n; i++) { f3[i] = h*f3[i]; } } /*************************************************************************** *********/ void F4(double t, const double *y, double u, double h, int n, const double *f3, double *f4, double *arg2) { int i; for (i = 0; i < n; i++) { arg2[i] = y[i] + f3[i]; } deriv((t+h), arg2, u, f4); for (i = 0; i < n; i++) { f4[i] = h*f4[i]; } } /*************************************************************************** *********/ void rk_init(RK_STRUCT *rks, int n) { /* Allocate memory. */ if (!(rks->f1 = (double *) malloc(n*sizeof(double)))) { printf("Error allocating memory.\n"); exit(1); } if (!(rks->f2 = (double *) malloc(n*sizeof(double)))) { printf("Error allocating memory.\n"); exit(1); } if (!(rks->f3 = (double *) malloc(n*sizeof(double)))) { printf("Error allocating memory.\n"); exit(1); } if (!(rks->f4 = (double *) malloc(n*sizeof(double)))) { printf("Error allocating memory.\n"); exit(1); } if (!(rks->arg2 = (double *) malloc(n*sizeof(double)))) { printf("Error allocating memory.\n"); exit(1); } } /*************************************************************************** *********/ void rk_free(RK_STRUCT *rks) { /* Free allocated memory */ free(rks->f1); free(rks->f2); free(rks->f3); free(rks->f4); free(rks->arg2); } /*************************************************************************** *********/ #if 1 /* Undefine main() to call rk4.c from another program. */ /* Simulate the pendulum system with fuzzy controller. */ void main(void) { double u = 0, t = 0, h = 0.001, /* Integration step size. */ y[3] = {0.1,0.0,0.0}, /* Initial conditions. */ tmax = 3; /* Simulation time period. */ int n = 3; /* System order. */ RK_STRUCT rks; FILE *fp; if (!(fp = fopen("output","w"))) { printf("Can't open file.\n"); exit(1); } /* Initialization area. */ rk_init(&rks,n); fuzzy_init(&fuzzy_system); /* Simulate For tmax seconds. */ while(t>--->>CUT-HERE>>--------------------------------- %***************************************************** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % DP2.2(a) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Fuzzy Cruise Controller %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % By: Kevin Passino % Version: 2/6/97 % (uses some code fragments developed earlier % by B. Klinehoffer and K. Passino) % % This program is used to simulate the PI fuzzy controller for % the automobile cruise control problem in Design Problem 2.2(a). % The variable names are chosen accordingly. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear % Clear all variables in memory % Initialize vehicle parameters m=1300; % Mass of the vehicle Ar=0.3; % Aerodynamic drag tau=0.2; % Engine/brake time constant d=100; % Constant frictional force % Initialize parameters for the fuzzy controller nume=11; % Number of input membership functions for the e % universe of discourse numie=11; % Number of input membership functions for the integral % of e universe of discourse g1=1;,g2=.01;,g0=1000; % Scaling gains for tuning membership functions for % universes of discourse for e, integral of e (what we % sometimes call "int e" below), and u respectively % These can be tuned to try to improve the performance. we=0.2*(1/g1); % we is half the width of the triangular input membership % function bases (note that if you change g0, the base width % will correspondingly change so that we always end % up with uniformly distributed input membership functions) % Note that if you change nume you will need to adjust the % "0.2" factor if you want membership functions that % overlap in the same way. wie=0.2*(1/g2); % Similar to we but for the int e universe of discourse base=0.4*g0; % Base width of output membership fuctions of the fuzzy % controller % Place centers of membership functions of the fuzzy controller: % Centers of input membership functions for the e universe of % discourse of fuzzy controller (a vector of centers) ce=[-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1]*(1/g1); % Centers of input membership functions for the int e universe of % discourse of fuzzy controller (a vector of centers) cie=[-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1]*(1/g2); % This next matrix specifies the rules of the fuzzy controller. % The entries are the centers of the output membership functions. % This choice represents just one guess on how to synthesize % the fuzzy controller. Notice the regularity % of the pattern of rules (basiscally we are using a type of % saturated index adding). Notice that it is scaled by g0, the % output scaling factor, since it is a normalized rule-base. % The rule-base can be tuned to try to improve performance. rules=[-1 -1 -1 -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0; -1 -1 -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2; -1 -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4; -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6; -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8; -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1; -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1; -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1 1; -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1 1 1; -0.2 0 0.2 0.4 0.6 0.8 1 1 1 1 1; 0 0.2 0.4 0.6 0.8 1 1 1 1 1 1]*g0; % Next, we initialize the simulation: t=0; % Reset time to zero index=1; % This is time's index (not time, its index). tstop=30; % Stopping time for the simulation (in seconds) %tstop=60; % Stopping time for the simulation (in seconds), last part % of problem (a), i.e., test input 3 step=0.01; % Integration step size x=[18;197.2;20]; % Intial condition on state of the car %x=0*x; % Use a zero initial condition for the last part % of problem (a) (in conjunction with test input 3). % % Next, we start the simulation of the system. This is the main % loop for the simulation of fuzzy cruise control system. % while t <= tstop v(index)=x(1); % Output of the car (velocity) % Next, we define the reference input vd (desired speed as % specified in the problem). if t<=10, vd(index)=18; end % First, define "test input 1" if t>10, vd(index)=22; end %if t<=10, vd(index)=18; end % This is "test input 2" (ramp) %if t>10, vd(index)=vd(index-1)+(4/1500); end % Ramp up 4 in 15 sec. %if t>25, vd(index)=22; end %vd(index)=22; % Use this input for the last part of problem (a) % This is "test input 3." % Fuzzy controller calculations: % First, for the given fuzzy controller inputs we determine % the extent at which the error membership functions % of the fuzzy controller are on (this is the fuzzification part). ie_count=0;,e_count=0; % These are used to count the number of % non-zero mf certainities of e and int e e(index)=vd(index)-v(index); % Calculates the error input for the fuzzy controller b(index)=x(3); % Sets the value of the integral of e % The following if-then structure fills the vector mfe % with the certainty of each membership fucntion of e for the % current input e. We use triangular membership functions. if e(index)<=ce(1) % Takes care of saturation of the left-most % membership function mfe=[1 0 0 0 0 0 0 0 0 0 0]; % i.e., the only one on is the % left-most one e_count=e_count+1;,e_int=1; % One mf on, it is the % left-most one. elseif e(index)>=ce(nume) % Takes care of saturation % of the right-most mf mfe=[0 0 0 0 0 0 0 0 0 0 1]; e_count=e_count+1;,e_int=nume; % One mf on, it is the % right-most one else % In this case the input is on the middle part of the % universe of discourse for e % Next, we are going to cycle through the mfs to % find all that are on for i=1:nume if e(index)<=ce(i) mfe(i)=max([0 1+(e(index)-ce(i))/we]); % In this case the input is to the % left of the center ce(i) and we compute % the value of the mf centered at ce(i) % for this input e if mfe(i)~=0 % If the certainty is not equal to zero then say % that have one mf on by incrementing our count e_count=e_count+1; e_int=i; % This term holds the index last entry % with a non-zero term end else mfe(i)=max([0,1+(ce(i)-e(index))/we]); % In this case the input is to the % right of the center ce(i) if mfe(i)~=0 e_count=e_count+1; e_int=i; % This term holds the index of the % last entry with a non-zero term end end end end % The following if-then structure fills the vector mfie with the % certainty of each membership fucntion of the integral of e % for its current value (to understand this part of the code see the above % similar code for computing mfe). Clearly, it could be more efficient to % make a subroutine that performs these computations for each of % the fuzzy system inputs. if b(index)<=cie(1) % Takes care of saturation of left-most mf mfie=[1 0 0 0 0 0 0 0 0 0 0]; ie_count=ie_count+1; ie_int=1; elseif b(index)>=cie(numie) % Takes care of saturation of the right-most mf mfie=[0 0 0 0 0 0 0 0 0 0 1]; ie_count=ie_count+1; ie_int=numie; else for i=1:numie if b(index)<=cie(i) mfie(i)=max([0,1+(b(index)-cie(i))/wie]); if mfie(i)~=0 ie_count=ie_count+1; ie_int=i; % This term holds last entry % with a non-zero term end else mfie(i)=max([0,1+(cie(i)-b(index))/wie]); if mfie(i)~=0 ie_count=ie_count+1; ie_int=i; % This term holds last entry % with a non-zero term end end end end % The next loops calculate the crisp output using only the non- % zero premise of error,e, and integral of the error b. This cuts down % computation time since we will only compute % the contribution from the rules that are on (i.e., a maximum of % four rules for our case). Minimum is used for the premise % and implication. num=0; den=0; for k=(e_int-e_count+1):e_int % Scan over e indices of mfs that are on for l=(ie_int-ie_count+1):ie_int % Scan over int e indices of mfs that are on prem=min([mfe(k) mfie(l)]); % Value of premise membership function % This next calculation of num adds up the numerator for the % defuzzification formula. rules(k,l) is the output center % for the rule.base*(prem-(prem)^2/2 is the area of a symmetric % triangle with base width "base" and that is chopped off at % a height of prem (since we use minimum to represent the % implication). Computation of den is similar but without % rules(k,l). num=num+rules(k,l)*base*(prem-(prem)^2/2); den=den+base*(prem-(prem)^2/2); end end u(index)=num/den; % Crisp output of fuzzy controller that is the input % to the vehicle. % Next, the Runge-Kutta equations are used to find the next state. % Clearly, it would be better to use a Matlab "function" for % F (but here we do not so we can have only one program). % (This does not implement the exact equations in the book but % if you think about the problem a bit you will see that these % really should be accurate enough. Here, we are not calculating % several intermediate values of the controller output, only % one per integration step; we do this simply to save some % computations.) time(index)=t; F=[(1/m)*(-Ar*x(1)^2 - d + x(2)) ; (1/tau)*(-x(2)+u(index)) ; vd(index)-x(1) ]; k1=step*F; xnew=x+k1/2; F=[(1/m)*(-Ar*xnew(1)^2 - d + xnew(2)) ; (1/tau)*(-xnew(2)+u(index)) ; vd(index)-xnew(1) ]; k2=step*F; xnew=x+k2/2; F=[(1/m)*(-Ar*xnew(1)^2 - d + xnew(2)) ; (1/tau)*(-xnew(2)+u(index)) ; vd(index)-xnew(1) ]; k3=step*F; xnew=x+k3; F=[(1/m)*(-Ar*xnew(1)^2 - d + xnew(2)) ; (1/tau)*(-xnew(2)+u(index)) ; vd(index)-xnew(1) ]; k4=step*F; x=x+(1/6)*(k1+2*k2+2*k3+k4); % Calculated next state t=t+step; % Increments time index=index+1; % Increments the indexing term so that % index=1 corresponds to time t=0. end % This end statement goes with the first "while" statement % in the program % % Next, we provide plots of the input and output of the vehicle % along with the reference speed that we want to track. % It is also easy to plot the two inputs to the fuzzy controller, % the integral of e (b(t)) and e(t) if you would like since these % values are also saved by the program. % subplot(211) plot(time,v,'-',time,vd,'--') grid on xlabel('Time (sec)') title('Vehicle speed (solid) and desired speed (dashed)') subplot(212) plot(time,u) grid on xlabel('Time (sec)') title('Output of fuzzy controller (input to the vehicle)') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % End of program % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %***************************************************** >>--->>CUT-HERE>>--------------------------------- %***************************************************** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % DP2.2(b) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Fuzzy Cruise Controller %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % By: Kevin Passino % Version: 2/6/97 % (uses some code fragments developed earlier % by B. Klinehoffer and K. Passino) % % This program is used to simulate the PD fuzzy controller for % the automobile cruise control problem in Design Problem 2.2(b). % The variable names are chosen accordingly. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear % Clear all variables in memory % Initialize vehicle parameters m=1300; % Mass of the vehicle Ar=0.3; % Aerodynamic drag tau=0.2; % Engine/brake time constant d=100; % Constant frictional force % Initialize parameters for the fuzzy controller nume=11; % Number of input membership functions for the e % universe of discourse numc=11; % Number of input membership functions for the c % universe of discourse g1=10;,g2=.5;,g0=1000; % Scaling gains for tuning membership functions for % universes of discourse for e, change in e % (what we call c) and u respectively % These can be tuned to try to improve the performance. we=0.2*(1/g1); % we is half the width of the triangular input membership % function bases (note that if you change g0, the base width % will correspondingly change so that we always end % up with uniformly distributed input membership functions) % Note that if you change nume you will need to adjust the % "0.2" factor if you want membership functions that % overlap in the same way. wc=0.2*(1/g2); % Similar to we but for the c universe of discourse base=0.4*g0; % Base width of output membership fuctions of the fuzzy % controller % Place centers of membership functions of the fuzzy controller: % Centers of input membership functions for the e universe of % discourse of fuzzy controller (a vector of centers) ce=[-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1]*(1/g1); % Centers of input membership functions for the int e universe of % discourse of fuzzy controller (a vector of centers) cc=[-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1]*(1/g2); % This next matrix specifies the rules of the fuzzy controller. % The entries are the centers of the output membership functions. % This choice represents just one guess on how to synthesize % the fuzzy controller. Notice the regularity % of the pattern of rules (basiscally we are using a type of % saturated index adding). Notice that it is scaled by g0, the % output scaling factor, since it is a normalized rule-base. % The rule-base can be tuned to try to improve performance. rules=[-1 -1 -1 -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0; -1 -1 -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2; -1 -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4; -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6; -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8; -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1; -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1; -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1 1; -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1 1 1; -0.2 0 0.2 0.4 0.6 0.8 1 1 1 1 1; 0 0.2 0.4 0.6 0.8 1 1 1 1 1 1]*g0; % Next, we initialize the simulation: t=0; % Reset time to zero index=1; % This is time's index (not time, its index). tstop=30; % Stopping time for the simulation (in seconds) %tstop=60; % Stopping time for the simulation (in seconds), last part % of problem (a), i.e., test input 3 step=0.01; % Integration step size eold=0; % Initialize the past value of the error (for use % in computing the derivative of the error e) x=[18;197.2]; % Intial condition on state of the car %x=0*x; % Use a zero initial condition for the last part % of problem (a) (in conjunction with test input 3). % % Next, we start the simulation of the system. This is the main % loop for the simulation of fuzzy cruise control system. % while t <= tstop v(index)=x(1); % Output of the car (velocity) % Next, we define the reference input vd (desired speed as % specified in the problem). if t<=10, vd(index)=18; end % First, define "test input 1" if t>10, vd(index)=22; end %if t<=10, vd(index)=18; end % This is "test input 2" (ramp) %if t>10, vd(index)=vd(index-1)+(4/1500); end % Ramp up 4 in 15 sec. %if t>25, vd(index)=22; end %vd(index)=22; % Use this input for the last part of problem (a) % This is "test input 3." % Fuzzy controller calculations: % First, for the given fuzzy controller inputs we determine % the extent at which the error membership functions % of the fuzzy controller are on (this is the fuzzification part). c_count=0;,e_count=0; % These are used to count the number of % non-zero mf certainities of e and c e(index)=vd(index)-v(index); % Calculates the error input for the fuzzy controller c(index)=(e(index)-eold)/step; % Sets the value of the c eold=e(index); % Save the past value of e for use in the above % computation the next time around the loop % The following if-then structure fills the vector mfe % with the certainty of each membership fucntion of e for the % current input e. We use triangular membership functions. if e(index)<=ce(1) % Takes care of saturation of the left-most % membership function mfe=[1 0 0 0 0 0 0 0 0 0 0]; % i.e., the only one on is the % left-most one e_count=e_count+1;,e_int=1; % One mf on, it is the % left-most one. elseif e(index)>=ce(nume) % Takes care of saturation % of the right-most mf mfe=[0 0 0 0 0 0 0 0 0 0 1]; e_count=e_count+1;,e_int=nume; % One mf on, it is the % right-most one else % In this case the input is on the middle part of the % universe of discourse for e % Next, we are going to cycle through the mfs to % find all that are on for i=1:nume if e(index)<=ce(i) mfe(i)=max([0 1+(e(index)-ce(i))/we]); % In this case the input is to the % left of the center ce(i) and we compute % the value of the mf centered at ce(i) % for this input e if mfe(i)~=0 % If the certainty is not equal to zero then say % that have one mf on by incrementing our count e_count=e_count+1; e_int=i; % This term holds the index last entry % with a non-zero term end else mfe(i)=max([0,1+(ce(i)-e(index))/we]); % In this case the input is to the % right of the center ce(i) if mfe(i)~=0 e_count=e_count+1; e_int=i; % This term holds the index of the % last entry with a non-zero term end end end end % The following if-then structure fills the vector mfc with the % certainty of each membership fucntion of the c % for its current value (to understand this part of the code see the above % similar code for computing mfe). Clearly, it could be more efficient to % make a subroutine that performs these computations for each of % the fuzzy system inputs. if c(index)<=cc(1) % Takes care of saturation of left-most mf mfc=[1 0 0 0 0 0 0 0 0 0 0]; c_count=c_count+1; c_int=1; elseif c(index)>=cc(numc) % Takes care of saturation of the right-most mf mfc=[0 0 0 0 0 0 0 0 0 0 1]; c_count=c_count+1; c_int=numc; else for i=1:numc if c(index)<=cc(i) mfc(i)=max([0,1+(c(index)-cc(i))/wc]); if mfc(i)~=0 c_count=c_count+1; c_int=i; % This term holds last entry % with a non-zero term end else mfc(i)=max([0,1+(cc(i)-c(index))/wc]); if mfc(i)~=0 c_count=c_count+1; c_int=i; % This term holds last entry % with a non-zero term end end end end % The next loops calculate the crisp output using only the non- % zero premise of error,e, and c. This cuts down % computation time since we will only compute % the contribution from the rules that are on (i.e., a maximum of % four rules for our case). Minimum is used for the premise % and implication. num=0; den=0; for k=(e_int-e_count+1):e_int % Scan over e indices of mfs that are on for l=(c_int-c_count+1):c_int % Scan over c indices of mfs that are on prem=min([mfe(k) mfc(l)]); % Value of premise membership function % This next calculation of num adds up the numerator for the % defuzzification formula. rules(k,l) is the output center % for the rule.base*(prem-(prem)^2/2 is the area of a symmetric % triangle with base width "base" and that is chopped off at % a height of prem (since we use minimum to represent the % implication). Computation of den is similar but without % rules(k,l). num=num+rules(k,l)*base*(prem-(prem)^2/2); den=den+base*(prem-(prem)^2/2); end end u(index)=num/den; % Crisp output of fuzzy controller that is the input % to the vehicle. % Next, the Runge-Kutta equations are used to find the next state. % Clearly, it would be better to use a Matlab "function" for % F (but here we do not so we can have only one program). % (This does not implement the exact equations in the book but % if you think about the problem a bit you will see that these % really should be accurate enough. Here, we are not calculating % several intermediate values of the controller output, only % one per integration step; we do this simply to save some % computations.) time(index)=t; F=[(1/m)*(-Ar*x(1)^2 - d + x(2)) ; (1/tau)*(-x(2)+u(index)) ]; k1=step*F; xnew=x+k1/2; F=[(1/m)*(-Ar*xnew(1)^2 - d + xnew(2)) ; (1/tau)*(-xnew(2)+u(index)) ]; k2=step*F; xnew=x+k2/2; F=[(1/m)*(-Ar*xnew(1)^2 - d + xnew(2)) ; (1/tau)*(-xnew(2)+u(index)) ]; k3=step*F; xnew=x+k3; F=[(1/m)*(-Ar*xnew(1)^2 - d + xnew(2)) ; (1/tau)*(-xnew(2)+u(index)) ]; k4=step*F; x=x+(1/6)*(k1+2*k2+2*k3+k4); % Calculated next state t=t+step; % Increments time index=index+1; % Increments the indexing term so that % index=1 corresponds to time t=0. end % This end statement goes with the first "while" statement % in the program % % Next, we provide plots of the input and output of the vehicle % along with the reference speed that we want to track. % It is also easy to plot the two inputs to the fuzzy controller, % the integral of e (b(t)) and e(t) if you would like since these % values are also saved by the program. % subplot(211) plot(time,v,'-',time,vd,'--') grid on xlabel('Time (sec)') title('Vehicle speed (solid) and desired speed (dashed)') subplot(212) plot(time,u) grid on xlabel('Time (sec)') title('Output of fuzzy controller (input to the vehicle)') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % End of program % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %***************************************************** >>--->>CUT-HERE>>--------------------------------- %***************************************************** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % E5.1(a) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % This program implements the batch least squares example % for data set G that is used in Chapter 5 of the book. % Kevin Passino, 1/14/97 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Initialize variables clear % The data set G M=3; X=[0 2 3; 2 4 6]; % In X the columns are the input vectors and the outputs are Y=[1; 5; 6]; n=2; R=2; % Next pick the input membership function centers to % lie in the middle of the training data at some type of % regular spacing - this is just a guess at where they should % be. C=[1.5 3; 3 5]; % Another set of possible input membership function centers % is given next. For this one we simply pick the centers to % be at the first R training points - again, this is just a % guess at where they should be and for some training data sets % if the first R training data pairs are not spread out properly % it may not be a good choice. For now, this option is commented % out - you can test it if you would like. You get similar % results to what we get for the choice of C above - and % hence draw similar conclusions. % %C=[0 2; % 2 4]; % Just pick one spread for the membership functions sigma=2; % Compute the xi (chi) vectors for the given data pairs for l=1:M % for all the training data for i=1:R mu(i,l)=1; % Initialization for j=1:n mu(i,l)=mu(i,l)*exp(-.5*((X(j,l)-C(j,i))/sigma)^2); % Numerator of xi % for lth training data pair for the % ith rule end end den(l)=sum(mu(:,l)); % Denominator of xi for lth training data xi(:,l)=mu(:,l)/den(l); % Compute xi vector for lth training data pair end % The above loop returns xi where each column is % the greek xi from the book for the lth data pair. Hence, Phi=xi'; % Set the weighting matrix in case you want to use weighted batch % least squares. W=eye(M); % This choice is for the non-weighted case. %W=diag([10 1 1]); % This weights the matching of the first data pair % more important. % Next, compute the least squares solution for the output membership % function centers. thetahat=inv(Phi'*W*Phi)*Phi'*W*Y % Test on the fuzzy system (at training data pairs and three other % test points). Xhat=[X [1;2] [2.5; 5] [4; 7]]; MGamma=6; % MGamma is the number of test points for l=1:MGamma % Next compute the fuzzy system - use the same approach as above % by computing the xi vector then multiplying it by the computed % parameter vector to form the sum in the numerator. for i=1:R muf(i,l)=1; % Initialization for j=1:n muf(i,l)=muf(i,l)*exp(-.5*((Xhat(j,l)-C(j,i))/sigma)^2); % Numerator of xi % for input test point end end denf(l)=sum(muf(:,l)); % Denominator of xi for lth training data xif(:,l)=muf(:,l)/denf(l); output(l)=thetahat'*xif(:,l); % Compute the output for the test point end % Next, display the input and output Xhat output %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %***************************************************** >>--->>CUT-HERE>>--------------------------------- %***************************************************** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % E5.1(c) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % This program implements the recursive least squares example % for data set G that is used in Chapter 5 of the book. % Kevin Passino, 1/14/97 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Initialize variables clear % The data set G M=3; X=[0 2 3; 2 4 6]; Xsave=X; % Save this for testing the fuzzy system later % In X the columns are the input vectors and the outputs are Y=[1; 5; 6]; % For RLS we will cycle through the training data set G NRLS times. % Hence, we will create a new training data set for convenience. NRLS=20; tempx=X; tempy=Y; for k=2:NRLS X=[X tempx]; % Makes columns that are the input portions Y=[Y ; tempy]; % Makes one long column end % Set number of inputs and rules n=2; R=2; % Next pick the input membership function centers to % lie in the middle of the training data at some type of % regular spacing - this is just a guess at where they should % be. C=[1.5 3; 3 5]; % Another set of possible input membership function centers % is given next. For this one we simply pick the centers to % be at the first R training points - again, this is just a % guess at where they should be and for some training data sets % if the first R training data pairs are not spread out properly % it may not be a good choice. For now, this option is commented % out - you can test it if you would like. You get similar % results to what we get for the choice of C above - and % hence draw similar conclusions. % %C=[0 2; % 2 4]; % Just pick one spread for the membership functions sigma=2; % Next set the initial condition for the RLS thetahatpast=[2; 5.5]; % This is thetahat at time zero %thetahatpast=[.3646; 8.1779]; % Another initialization option - the % result from BLS thetahatinitial=thetahatpast; % Save it for later use alpha=2000; Ppast=alpha*eye(R); % This is P at time zero lambda=1; % Use non-weighted least squares (but can change this % to, e.g., lambda=0.9 to add a forgetting factor. % Next, we have the main loop where we use RLS to train % the fuzzy system. for k=1:NRLS*M % For all the training data % Compute the xi (chi) vectors for the given data pairs. for i=1:R mu(i,k)=1; % Initialization for j=1:n mu(i,k)=mu(i,k)*exp(-.5*((X(j,k)-C(j,i))/sigma)^2); % Numerator of xi % for kth training data pair for the % ith rule end end den(k)=sum(mu(:,k)); % Denominator of xi for kth training data xi(:,k)=mu(:,k)/den(k); % Compute xi vector for kth training data pair % Next, we specify the RLS equations % P is P(k), i.e., P at time k. Ppast is P(k-1). P=(1/lambda)*(eye(R)-Ppast*xi(:,k)*inv(lambda + ... xi(:,k)'*Ppast*xi(:,k))*xi(:,k)')*Ppast; thetahat(:,k)=thetahatpast+... P*xi(:,k)*(Y(k)-xi(:,k)'*thetahatpast); Ppast=P; thetahatpast=thetahat(:,k); end thetahat=[thetahatinitial thetahat]; % Save whole vector for plotting k=0:NRLS*M; % Set time vector for plotting % Plot the RLS estimates vs. time plot(k,thetahat) xlabel('Time step, k') title('RLS: thetahat_1 and thetahat_2') % Next, we let the RLS solution for the output membership % function centers be the last computed value of thetahat. % Print out the last estimated value and the corresponding P % matrix. thetahatresult=thetahatpast Ppast % Test the fuzzy system at training data pairs and three other % test points. Xhat=[Xsave [1;2] [2.5; 5] [4; 7]]; MGamma=6; % MGamma is the number of test points for l=1:MGamma % Next compute the fuzzy system - use the same approach as above % by computing the xi vector then multiplying it by the computed % parameter vector to form the sum in the numerator. for i=1:R muf(i,l)=1; % Initialization for j=1:n muf(i,l)=muf(i,l)*exp(-.5*((Xhat(j,l)-C(j,i))/sigma)^2); % Numerator of xi % for input test point end end denf(l)=sum(muf(:,l)); % Denominator of xi for lth input data xif(:,l)=muf(:,l)/denf(l); output(l)=thetahatresult'*xif(:,l); % Compute the output for % the test point end % Next, display the input and output Xhat output %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %***************************************************** >>--->>CUT-HERE>>--------------------------------- %***************************************************** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % E5.1(e) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % This program implements the gradient method example % for data set G that is used in Chapter 5 of the book. % Kevin Passino, 1/15/97 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Initialize variables clear % The data set G M=3; X=[0 2 3; 2 4 6]; Xsave=X; % Save this for testing the fuzzy system later % In X the columns are the input vectors and the outputs are Y=[1; 5; 6]; % For the gradient method we will cycle through the training data % set G Ngrad times. Hence, we will create a new training data % set for convenience. Ngrad=40; tempx=X; tempy=Y; for k=2:Ngrad X=[X tempx]; % Makes columns that are the input portions Y=[Y ; tempy]; % Makes one long column end % Set number of inputs and rules n=2; R=2; % Next set the initial conditions for the Gradient method % First, pick the *initial* values of the consequent linear (affine) % functions. These are in the A matrix with the first n columns % the terms that multiply x^1 to x^n and the last column is the % bias term (the zeroth term). Each row of A corresponds to a % different rule. Now, the initialization can be done in many ways. % You could use batch least squares to pick them for your initial guesses % for the input membership function centers and spreads below. You % could have some guesses based on insights from the problem. Here, % we simply take a guess for these parameters and hope the % algorithm converges. Apast=[1 1 1; 1 1 1]; % Next, pick the *initial* input membership function centers to % lie in the middle of the training data at some type of % regular spacing - this is just a guess at where they should % be. Cpast=[1.5 3; 3 5]; % Another set of possible initial input membership function centers % is given next. For this one we simply pick the centers to % be at the first R training points - again, this is just a % guess at where they should be and for some training data sets % if the first R training data pairs are not spread out properly % it may not be a good choice. For now, this option is commented % out - you can test it if you would like. You get similar % results to what we get for the choice of C above - and % hence draw similar conclusions. % %Cpast=[0 2; % 2 4]; % Finally, pick an initial spread for all the input membership functions. % Here, we simply pick all the intial spreads to be the same. Sigmapast=[2 2; 2 2]; % Note that these "past" values correspond to time zero, but this will % be indexed with time one since Matlab can't use the zero index. % Next, we choose the step sizes for the algorithm. lambda_2=1; % Step size for the c parameters in C lambda_3=1; % Step size for the sigma parameters in Sigma lambda_4=0.01; % Step size for the a_i,j parameters in A % Next, we have the main loop where we use RLS to train % the fuzzy system. for k=1:Ngrad*M % For all the training data % Compute the xi (chi) vectors for the given data pair. for i=1:R mu(i,k)=1; % Initialization g(i,k)=Apast(i,n+1); % Stored the bias term for the consequent as % the last element of the rows of G for j=1:n mu(i,k)=mu(i,k)*exp(-.5*((X(j,k)-Cpast(j,i))/Sigmapast(j,i))^2); % Numerator of xi % for kth training data pair for the % ith rule g(i,k)=g(i,k)+Apast(i,j)*X(j,k); % This sums up the consequent end end den(k)=sum(mu(:,k)); % Denominator of xi for kth training data xi(:,k)=mu(:,k)/den(k); % Compute xi vector for kth training data pair f(k)=g(:,k)'*xi(:,k); % Compute the output for % for the current training data epsilon(k)=f(k)-Y(k); % For use in the gradient updates e(k)=0.5*(epsilon(k))^2; % Compute the approximation error % While often you would, we do not use % this to help specify a termination % condition. Here, we simply train % for a fixed number of steps. % Next, we specify the Gradient equations % First, set a value for which no input membership function spread can % decrease below - call this value sigmabar sigmabar=.01; for i=1:R for j=1:n A(i,j)=Apast(i,j)-lambda_4*epsilon(k)*xi(i,k)*X(j,k); % Update consequent parameters C(j,i)=Cpast(j,i)-lambda_2*epsilon(k)*((g(i,k)-f(k))/den(k))*... mu(i,k)*((X(j,k)-Cpast(j,i))/(Sigmapast(j,i))^2); % Update input mf centers Sigma(j,i)=Sigmapast(j,i)-lambda_3*epsilon(k)*((g(i,k)-f(k))/den(k))*... mu(i,k)*((X(j,k)-Cpast(j,i))^2/(Sigmapast(j,i))^3); % Update input mf spreads if Sigma(j,i)>--->>CUT-HERE>>--------------------------------- %***************************************************** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % E5.1(f) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This program performs the necessary computations % for the Chapter 5 example for clustering with % optimal output predefuzzification. % K. Passino, 12/19/96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Initialize variables clear % The data set G M=3; X=[0 2 3; 2 4 6]; % In X the columns are the input vectors Y=[1; 5; 6]; % Initial cluster centers R=2; V=[1.89 2.47; 3.76 4.76]; % Here, the columns of V are the cluster centers % Set other paramters for the clustering algorithm e_c=0.001; % This is epsilon_c in the algorithm m=2; % This is the fuzziness factor num_iter=0; % Counter for number of iterations flg=0; % Termination flag used in the program % First loop is for then number of iterations % in the clustering method (flg=1 when meet % the termination criterion at the end of the % program). while flg==0 num_iter=num_iter+1; %Count the number of iterations % Compute the mu_ij values for i=1:M for j=1:R z=0; for k=1:R z=z+(((X(:,i)-V(:,j))'*(X(:,i)-V(:,j)))/((X(:,i)-V(:,k))'*(X(:,i)-V(:,k))))^(1/(m-1)); end mu(i,j)=1/z; end end mu % Compute the new cluster centers for j=1:R num=[0 ; 0]; den=0; for i=1:M num=num+X(:,i)*(mu(i,j)^m); den=den+(mu(i,j)^m); end Vnew(:,j)=num/den; end Vnew % Check for how much the cluster centers moved for i=1:R test(i)=sqrt((V(:,i)-Vnew(:,i))'*(V(:,i)-Vnew(:,i))); end test % Use this to test if the cluster centers have both moved % less than e_c and if they have moved very little % then terminate if test(1)<=e_c & test(2)<=e_c, flg=1; end % Before repeating the loop make the newly computed % cluster centers the old ones V=Vnew; end % Print out the number of iterations num_iter % Next, we perform least squares computations for determining % the parameters of the linear (affine) consequent functions. Xhat=[ 1 1 1; X]; % Dhat is the vector of elements % that we will place on the diagonal for j=1:R Dhat(:,j)=[mu(1,j) mu(2,j) mu(3,j)]'; end % Next we compute the weighted least squares solution for j=1:R A(:,j)=inv(Xhat*(diag(Dhat(:,j))^2)*Xhat')*Xhat*(diag(Dhat(:,j))^2)*Y; end A % Next we test the Takagi-Sugeno fuzzy system with % an input point [1 2]' x=[1; 2]; numf=0; denf=0; for j=1:R z=0; for k=1:R z=z+(((x-V(:,j))'*(x-V(:,j)))/((x-V(:,k))'*(x-V(:,k))))^(1/(m-1)); end muh(j)=1/z; numf=numf+(A(1,j) + A(2,j)*x(1) + A(3,j)*x(2))*muh(j); denf=denf+muh(j); end % Next, display the output output=numf/denf %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %***************************************************** >>--->>CUT-HERE>>--------------------------------- %***************************************************** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % E5.1(f) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This program performs the necessary computations % for the Chapter 5 example for clustering with % optimal output predefuzzification. This is a modified % program that tests having 4 i-o data pairs. % K. Passino, 1/7/97 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Initialize variables clear % The data set G M=4; X=[0 2 3 4; 2 4 6 5]; % In X the columns are the input vectors Y=[1; 5; 6; 5.5]; % Initial cluster centers R=2; V=[1.89 2.47; 3.76 4.76]; % Here, the columns of V are the cluster centers % Set other paramters for the clustering algorithm e_c=0.001; % This is epsilon_c in the algorithm m=2; % This is the fuzziness factor num_iter=0; % Counter for number of iterations flg=0; % Termination flag used in the program % First loop is for then number of iterations % in the clustering method (flg=1 when meet % the termination criterion at the end of the % program). while flg==0 num_iter=num_iter+1; %Count the number of iterations % Compute the mu_ij values for i=1:M for j=1:R z=0; for k=1:R z=z+(((X(:,i)-V(:,j))'*(X(:,i)-V(:,j)))/((X(:,i)-V(:,k))'*(X(:,i)-V(:,k))))^(1/(m-1)); end mu(i,j)=1/z; end end mu % Compute the new cluster centers for j=1:R num=[0 ; 0]; den=0; for i=1:M num=num+X(:,i)*(mu(i,j)^m); den=den+(mu(i,j)^m); end Vnew(:,j)=num/den; end Vnew % Check for how much the cluster centers moved for i=1:R test(i)=sqrt((V(:,i)-Vnew(:,i))'*(V(:,i)-Vnew(:,i))); end test % Use this to test if the cluster centers have both moved % less than e_c and if they have moved very little % then terminate if test(1)<=e_c & test(2)<=e_c, flg=1; end % Before repeating the loop make the newly computed % cluster centers the old ones V=Vnew; end % Print out the number of iterations num_iter % Next, we perform least squares computations for determining % the parameters of the linear (affine) consequent functions. Xhat=[ 1 1 1 1; X]; % Dhat is the vector of elements % that we will place on the diagonal for j=1:R Dhat(:,j)=[mu(1,j) mu(2,j) mu(3,j) mu(4,j)]'; end % Next we compute the weighted least squares solution for j=1:R A(:,j)=inv(Xhat*(diag(Dhat(:,j))^2)*Xhat')*Xhat*(diag(Dhat(:,j))^2)*Y; end A % Next we test the Takagi-Sugeno fuzzy system with % an input point [1 2]' x=[1; 2]; numf=0; denf=0; for j=1:R z=0; for k=1:R z=z+(((x-V(:,j))'*(x-V(:,j)))/((x-V(:,k))'*(x-V(:,k))))^(1/(m-1)); end muh(j)=1/z; numf=numf+(A(1,j) + A(2,j)*x(1) + A(3,j)*x(2))*muh(j); denf=denf+muh(j); end % Next, display the output output=numf/denf %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %***************************************************** >>--->>CUT-HERE>>--------------------------------- %***************************************************** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % E6.1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Fuzzy Model Reference Learning Control System Simulator %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % By: Brian Klinehoffer and Kevin Passino % Version: 7/6/98 % In earlier version there was a typo in the part of the code % that computed mfye where "e" should have been "ye". % % This program is used to simulate the "fuzzy model reference % learning controller" (FMRLC) for a simple plant. The objective % here is to illustrate the basic issues in how to simulate this % adaptive fuzzy controller and to show how it will perform under % various conditions. % % The plant is a second order system in the form of % k_p/(s^2 + 2*zeta_p*w_p*s + w_p^2) where k_p is the gain in the % numerator, zeta_p is the damping ratio and w_p is the undamped % natural frequency. We encourage the user to change these parameters % (e.g., even on-line) to see how the FMRLC can adapt to the % changes in the plant that such a variation represents. We use % the simple linear system for the sake of illustration. We have % found the FMRLC particularly well-suited for a variety of nonlinear % plants and it is for these that the FMRLC shows some distinct % advantages over conventional methods (for some applications). % Suppose that we name the input to the plant u and its output y. % % The goal of the FMRLC is to make the closed-loop system behave like % a given "reference model." Here, our reference model is a % first order system that represents how we would like the system % to behave. For instance, the plant may have unknown pole positions % with a low zeta_p value so that the plant is oscillatory. The % reference model could specify that we would like the pole of % the closed-loop system to be further out into the left-half plane % and have a higher zeta so that it will be faster. The input to % the reference model is the reference input r and its output is ym. % % There are four components to the FMRLC: the plant, reference model, % fuzzy controller that will be tuned, and "fuzzy inverse % model." We will not fully explain the rationale for the choice of % the parameters for each of these in this program. Our objective % here is simply to show how to program the FMRLC once you have a % basic understanding of its operation. % % The fuzzy controller here has two inputs and one output. One input % is the error e=r-y and the other is the change in error % c=(e(kT)-e(kt-T))/T where T=step (in our program). The fuzzy % controller uses 11 uniformly distributed triangular membership % functions for each of its input universes of discourse. We use % "linguistic-numeric" indices of 1,2,3,...,11 to name the membership % functions (similarly for the fuzzy inverse model). % The fuzzy controller uses minimum to represent the premise % and implication, and COG defuzzification. The output membership % function centers of the fuzzy controller are what is tuned by the % learning mechanism of the FMRLC. Clearly, all these choices for % the structure and parameters are simply for illustrative purposes; % other choices may be more appropriate for other plants and control % objectives. % % The learning mechanism consists of two parts: the fuzzy inverse % model and the knowledge-base modifier. The fuzzy inverse model % is simply a fuzzy controller that operates in the adaptation loop. % Its two inputs are ye=ym - y and yc=(ye(kT)-ye(kT-T))/T (again % T=step). It uses 11 uniformly distributed triangular membership % functions for each of its input universes of discourse. It uses % triangular shaped output membership functions. It uses minimum % to represent the premise and implication, and COG defuzzification. % The output of the fuzzy inverse model is p and it indicates how % much to move the centers of the output membership functions. % The knowledge-base modifier is the mechanism that actually moves % the fuzzy controller output membership function centers. Here, % we use the fact that either 1, 2, or 4 rules of the fuzzy % controller will be on at any one time. Our update formula is % c_i(kT)=c_i(kT-dT)+p(kT) where c_i is the output center for % the ith rule and the index i is selected from the set of rules % that were on d steps before (at most 4 rules) (we use different % notation in the program below). In this way the rules that % contributed to the current error y_e are the ones that % are modified by the learning mechanism. It is important to note % that while for the simple linear system we do not design a very % innovative fuzzy inverse model, real-world applications will % dictate the need for much more careful tuning/design of the % rule-base and membership functions for the inverse model. The % one we have chosen works adequately for this simple problem. % Next, note that in the code we have taken certain shortcuts % soley for illustrative purposes. For instance, the code for the % fuzzy controller that is tuned is very similar to that of the % fuzzy inverse model; while this works for this simple application % clearly it will not always be possible to use such similar % fuzzy systems (e.g., ones that have the same number of inputs, % rules, and membership functions). This completes the description % of the FMRLC. % % The FMRLC is simulated with a 4th order Runge-Kutta technique. One % must be careful in the choice of the integration step size. For % instance, if you change the a_p parameter to be very large, which % corresponds to a fast system, then you will need to make the value % of the integration step size smaller so that the simulation is % accurate. This completes our brief description of how this program % operates; however, comments are included throughout the program to % further clarify its operation. If you want more details on how the % FMRLC works see the following references: % % Layne J.R., Passino K.M., "Fuzzy Model Reference Learning Control % for Cargo Ship Steering", IEEE Control Systems Magazine, Vol. 13, % No. 6, pp. 23-34, Dec. 1993. % % Layne J.R., Passino K.M., "Fuzzy Model Reference Learning Control," % To appear in the Journal of Intelligent and Fuzzy Systems, 1996. % % Layne J.R., Passino K.M., Yurkovich S., "Fuzzy Learning Control % for Anti-Skid Braking Systems", IEEE Trans. on Control Systems % Technology, Vol. 1, No. 2, pp. 122-129, June 1993. % % Kwong W.A., Passino K.M., Lauknonen E.G., Yurkovich S., ÒExpert % Supervision of Fuzzy Learning Systems for Fault Tolerant Aircraft % ControlÓ, Special Issue on Fuzzy Logic in Engineering Applications, % Proceedings of the Inst. of Electrical and Electronics Engineers, % Vol. 83, No.3, pp. 466-483, March 1995. % % Moudgal V.G., Kwong W.A., Passino K.M., and Yurkovich S., ÒFuzzy % Learning Control for a Flexible-Link RobotÓ, IEEE Transactions on % Fuzzy Systems, Vol. 3, No. 2, pp. 199-210, May 1995. % % Kwong W.A., Passino K.M., "Dynamically Focused Fuzzy Learning % Control", IEEE Trans. on Systems, Man, and Cybernetics, Vol. 26, % No. 1, pp. 53-74, Feb. 1996. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear % Clear all variables in memory % Initialize variables eold=0; % Intial condition used to calculate c rold=0; % Intial condition used to calculate r yeold=0; % Intial condition used to calculate yc ymold=0; % Initial condition for the first order reference model % Next, initialize parameters for the fuzzy controller nume=11; % Number of input membership functions for the e % universe of discourse numc=11; % Number of input membership functions for the c % universe of discourse ge=1/2;,gc=1/2;,gu=5; % Scaling gains for tuning membership functions for % universes of discourse for e, c and u respectively % These are tuned to improve the performance of the FMRLC we=0.2*(1/ge); % we is half the width of the triangular input membership % function bases (note that if you change ge, the base width % will correspondingly change so that we always end % up with uniformly distributed input membership functions) % Note that if you change nume you will need to adjust the % "0.2" factor if you want membership functions that % overlap in the same way. wc=0.2*(1/gc); % Similar to we but for the c universe of discourse base=0.4*gu; % Base width of output membership fuctions of the fuzzy % controller % Place centers of membership functions of the fuzzy controller: % Centers of input membership functions for the e universe of % discourse for of fuzzy controller (a vector of centers) ce=[-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1]*(1/ge); % Centers of input membership functions for the c universe of % discourse for of fuzzy controller (a vector of centers) cc=[-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1]*(1/gc); % This next matrix determines the rules of the fuzzy controller. % The entries are the centers of the output membership functions. % This choice represents our "best guess" on how to synthesize % the fuzzy controller for the nominal system (with a little % work you could certainly do better). Notice the regularity % of the pattern of rules (basiscally we are using a type of % saturated index adding). Notice that it is scaled by gu, the % output scaling factor, since it is a normalized rule-base % % Also, there is the parameter gf that is either =0 or =1. It % mulitplies the rule-base and is =1 if you want the FMRLC to be % initialized with a decent guess at what the fuzzy controller % that it is going to tune should be. If gf=0 then the rule-base % has 121 rules with output membership functions that are triangular % with base widths base=0.4*gu and centers at zero. Hence, % if gf=0 this represents that the fuzzy controller does not % know very much at all about how to control the plant initially. % As the FMRLC learns it will fill in the rule-base to % try to achieve good behavior. gf=0; fuzzyrules=[-1 -1 -1 -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0; -1 -1 -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2; -1 -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4; -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6; -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8; -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1; -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1; -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1 1; -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1 1 1; -0.2 0 0.2 0.4 0.6 0.8 1 1 1 1 1; 0 0.2 0.4 0.6 0.8 1 1 1 1 1 1]*gu*gf; % Next, we define some parameters for the fuzzy inverse model gye=1/2;,gyc=1/2; % Scaling gains for the error and change in error for % the inverse model % These are tuned to improve the performance of the FMRLC gp=0.2; % Scaling gain for the output of inverse model. If you let gp=0 then % you turn off the learning mechanism and hence the FMRLC reduces to % a direct fuzzy controller that is not tuned. Hence, from this % program you can also see how to code a non-adaptive fuzzy controller. % A reasonable value for it is gp=.2 but you should try other values. numye=11; % Number of input membership functions for the ye % universe of discourse numyc=11; % Number of input membership functions for the yc % universe of discourse wye=0.2*(1/gye); % Sets the width of the membership functions for % ye from center to extremes wyc=0.2*(1/gyc); % Sets the width of the membership functions for % yc from center to extremes invbase=0.4*gp; % Sets the base of the output membership functions % for the inverse model % Place centers of inverse model membership functions % For error input for learning Mechanism cye=[-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1]*(1/gye); % For change in error input for learning mechanism cyc=[-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1]*(1/gyc); % The next matrix contains the rule-base matrix for the fuzzy % inverse model. Notice that for simplicity we choose it to have % the same structure as the rule-base for the fuzzy controller. % While this will work for the control of the simple first order % linear system for many nonlinear systems a different structure % will be needed for the rule-base. Again, the entries are % the centers of the output membership functions, but now for % the fuzzy inverse model. inverrules=[-1 -1 -1 -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0; -1 -1 -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2; -1 -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4; -1 -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6; -1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8; -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1; -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1; -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1 1; -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1 1 1; -0.2 0 0.2 0.4 0.6 0.8 1 1 1 1 1; 0 0.2 0.4 0.6 0.8 1 1 1 1 1 1]*gp; % Next, we set up some parameters/variables for the % knowledge-base modifier d=1; % This sets the number of steps the knowledge-base modifier looks % back in time. For this program it must be an integer % less than or equal to 10 (but this is easy to make larger) % The next four vectors are used to store the information about % which rules were on 1 step in the past, 2 steps in the past, ...., % 10 steps in the past (so that picking 0<= d <= 10 can be used). meme_int=[0 0 0 0 0 0 0 0 0 0]; % sets up the vector to store up to 10 values of e_int meme_count=[0 0 0 0 0 0 0 0 0 0]; % sets up the vector to store up to 10 values of e_count memc_int=[0 0 0 0 0 0 0 0 0 0]; % sets up the vector to store up to 10 values of c_int memc_count=[0 0 0 0 0 0 0 0 0 0]; % sets up the vector to store up to 10 values of c_count % % Next, we intialize the simulation of the closed-loop system. % k_p=1; % The numerator of the plant. Change this value to study % the ability of the FMRLC to control other plants. Also, % you can make this a time-varying parameter. zeta_p=.707; % Damping ratio for the second order plant (could change this % to see how the system will adapt to it) w_p=1; % Undamped natural frequency for the plant (could change this % to see how the system will adapt to it) k_r=1; % The numerator of the reference model. Change this value to study % the ability of the FMRLC to meet other performance specifications. a_r=1; % The value of -a_r is the pole position for the reference model. % Change this value to study the ability of the FMRLC to meet % other performance specifications (e.g., a faster response). t=0; % Reset time to zero index=1; % This is time's index (not time, its index). tstop=64; % Stopping time for the simulation (in seconds) step=0.01; % Integration step size x=[0;0]; % Intial condition on state of the plant % Need a state space representation for the plant. Since our % plant is linear we use the standard form of xdot=Ax+Bu, y=Cx+Du % Matrix A of state space representation of plant A=[0 1; -w_p^2 -2*zeta_p*w_p]; B=[0; 1]; % Matrix B of state space representation of plant C=[k_p 0]; % Matrix C of state space representation of plant % % Next, we start the simulation of the system. This is the main % loop for the simulation of the FMRLC. % while t <= tstop y(index)=C*x; % Output of the plant % Next, we define the reference input r as a sine wave r(index)=sin(.6*t); % The following sets up the reference model. Just for the sake of % illustration we use a discrete-time reference model. If you use % a continuous-time reference model you will have to augment % the state of the closed-loop system with the state(s) of the % reference model and hence update the state in the Runge-Kutta % equations above. % % For the reference model we use a first order transfer function % k_r/(s+a_r) but we use the bilinear transformation where we % replace s by (2/step)(z-1)/(z+1), then find the z-domain % representation of the reference model, then convert this % to a difference equation: ym(index)=(1/(2+a_r*step))*((2-a_r*step)*ymold+... k_r*step*(r(index)+rold)); ymold=ym(index); rold=r(index); % This saves the past value of the ym (r) so that we can use it % the next time around the loop % Now that we have simulated the next step for the plant and reference % model we will focus on the two fuzzy components. % First, for the given fuzzy controller inputs we determine % the extent at which the error membership functions % of the fuzzy controller are on (this is the fuzzification part). c_count=0;,e_count=0; % These are used to count the number of % non-zero mf certainities of e and c e=r(index)-y(index); % Calculates the error input for the fuzzy controller c=(e-eold)/step; % Calculates the change in error input for the fuzzy controller eold=e; % Saves the past value of e for use in the next time through the % loop % The following if-then structure fills the vector mfe % with the certainty of each membership fucntion of e for the % current input e if e<=ce(1) % Takes care of saturation of the left-most % membership function mfe=[1 0 0 0 0 0 0 0 0 0 0]; % i.e., the only one on is the % left-most one e_count=e_count+1;,e_int=1; % One mf on, it is the % left-most one. elseif e>=ce(nume) % Takes care of saturation % of the right-most mf mfe=[0 0 0 0 0 0 0 0 0 0 1]; e_count=e_count+1;,e_int=nume; % One mf on, it is the % right-most one else % In this case the input is on the middle part of the % universe of discourse for e % Next, we are going to cycle through the mfs to % find all that are on for i=1:nume if e<=ce(i) mfe(i)=max([0 1+(e-ce(i))/we]); % In this case the input is to the % left of the center ce(i) and we compute % the value of the mf centered at ce(i) % for this input e if mfe(i)~=0 % If the certainty is not equal to zero then say % that have one mf on by incrementing our count e_count=e_count+1; e_int=i; % This term holds the index last entry % with a non-zero term end else mfe(i)=max([0,1+(ce(i)-e)/we]); % In this case the input is to the % right of the center ce(i) if mfe(i)~=0 e_count=e_count+1; e_int=i; % This term holds the index of the % last entry with a non-zero term end end end end % Next we will save the number of mfs that are on and the pointer % e_int as to which rules were on. This vector of length % 10 saves the last 10 values of e_count and e_int as time % progresses (hence, it is a moving window). % These will be used by the FMRLC knowledge-base modifier. meme_count=[e_count meme_count(1:9)]; meme_int=[e_int meme_int(1:9)]; % The following if-then structure fills the vector mfc with the % certainty of each membership fucntion of c for the current % value of c (to understand this part of the code see the above % similar code for computing mfe) if c<=cc(1) % Takes care of saturation of left-most mf mfc=[1 0 0 0 0 0 0 0 0 0 0]; c_count=c_count+1; c_int=1; elseif c>=cc(numc) % Takes care of saturation of the right-most mf mfc=[0 0 0 0 0 0 0 0 0 0 1]; c_count=c_count+1; c_int=numc; else for i=1:numc if c<=cc(i) mfc(i)=max([0,1+(c-cc(i))/wc]); if mfc(i)~=0 c_count=c_count+1; c_int=i; % This term holds last entry % with a non-zero term end else mfc(i)=max([0,1+(cc(i)-c)/wc]); if mfc(i)~=0 c_count=c_count+1; c_int=i; % This term holds last entry % with a non-zero term end end end end % Next we will save the number of mfs that are on and the pointer % c_int as to which rules were on. This vector of length 10 % saves the last 10 values of c_count and e_int as time progresses % (hence, it is a moving window). These will be used by the FMRLC % knowledge-base modifier. memc_count=[c_count memc_count(1:9)]; memc_int=[c_int memc_int(1:9)]; % These for loops calculate the crisp output using only the non- % zero premise of error,e, and change in error,c. This cuts down % computation time since we will only compute % the contribution from the rules that are on (i.e., a maximum of % four rules for our case). Minimum is used for the premise % and implication. num=0; den=0; for k=(e_int-e_count+1):e_int % Scan over e indices of mfs that are on for l=(c_int-c_count+1):c_int % Scan over c indices of mfs that are on prem=min([mfe(k) mfc(l)]); % Value of premise membership function % This next calculation of num adds up the numerator for the % defuzzification formula. fuzzyrules(k,l) is the output center % for the rule.base*(prem-(prem)^2/2 is the area of a symmetric % triangle with base width "base" and that is chopped off at % a height of prem (since we use minimum to represent the % implication). Computation of den is similar but without % fuzzyrules(k,l). num=num+fuzzyrules(k,l)*base*(prem-(prem)^2/2); den=den+base*(prem-(prem)^2/2); end end u(index)=num/den; % Crisp output of fuzzy controller that is the input % to the plant % Next, we perform computations for the fuzzy inverse model. ye=ym(index)-y(index); % Calculates ye yc=(ye-yeold)/step; % Calculates yc yeold=ye; % Saves the value of ye for use the % next time ye_count=0;,yc_count=0; % Counts the non-zero certainties % of ye and yc similar to how we did % for the fuzzy controller % The following if-then structure fills the vector mfye with the % certainty of each membership fucntion of ye (similar to the % fuzzy controller). Notice that we use the same number of % input membership functions as in the fuzzy controller - % just for convenience - it does not have to be this way if ye<=cye(1) % Takes care of saturation mfye=[1 0 0 0 0 0 0 0 0 0 0]; ye_count=ye_count+1;,ye_int=1; elseif ye>=cye(numye) % Takes care of saturation mfye=[0 0 0 0 0 0 0 0 0 0 1]; ye_count=ye_count+1;,ye_int=numye; else for i=1:numye if ye<=cye(i) mfye(i)=max([0 1+(ye-cye(i))/wye]); if mfye(i)~=0 ye_count=ye_count+1; ye_int=i; % This term holds last entry with % a non-zero term end else mfye(i)=max([0,1+(cye(i)-ye)/wye]); if mfye(i)~=0 ye_count=ye_count+1; ye_int=i; % This term holds last entry with % a non-zero term end end end end % The following if-then structure fills the vector mfyc with the % certainty of each membership fucntion of yc if yc<=cyc(1) mfyc=[1 0 0 0 0 0 0 0 0 0 0]; yc_count=yc_count+1;,yc_int=1; elseif yc>=cyc(numyc) mfyc=[0 0 0 0 0 0 0 0 0 0 1]; yc_count=yc_count+1;,yc_int=numyc; else for i=1:numyc if yc<=cyc(i) mfyc(i)=max([0 1+(yc-cyc(i))/wyc]); if mfyc(i)~=0 yc_count=yc_count+1; yc_int=i; end else mfyc(i)=max([0,1+(cyc(i)-yc)/wyc]); if mfyc(i)~=0 yc_count=yc_count+1; yc_int=i; end end end end % These for loops calculate the crisp output of the inverse model % using only the non-zero premise of error,ye, and change in % error,yc. This cuts down computation time (similar to the % fuzzy controller). Minimum is used for the premise and the % implication. To understand this part of the code see the % similar code for the fuzzy controller. invnum=0; invden=0; for k=(ye_int-ye_count+1):ye_int for l=(yc_int-yc_count+1):yc_int prem=min([mfye(k) mfyc(l)]); invnum=invnum+inverrules(k,l)*invbase*(prem-(prem)^2/2); invden=invden+invbase*(prem-(prem)^2/2); end end % Next we compute the output of the fuzzy inverse model. % If you want to let gp=0 to test the fuzzy controller by itself % then you will have invden=0 and you will not be able to compute p. % To make it possible to let gp=0 we put in the following % if-then rule. if gp==0 p=0; else p=invnum/invden; % Crisp output of inverse model %if abs(p)<.05, p=0; end % robustification term end % Next we modify the centers of ouput membership functions % associated with the rules that were on d steps ago % The indexing sheme is similar to the one used in the % computation of the outputs of the fuzzy controller and % fuzzy inverse model. for k=(meme_int(d)-meme_count(d)+1):meme_int(d) for l=(memc_int(d)-memc_count(d)+1):memc_int(d) fuzzyrules(k,l)=fuzzyrules(k,l)+p; end end % Next, the Runge-Kutta equations to find the next state. time(index)=t; F=A*x+B*u(index); k1=step*F; xnew=x+k1/2; F=A*xnew+B*u(index); k2=step*F; xnew=x+k2/2; F=A*xnew+B*u(index); k3=step*F; xnew=x+k3; F=A*xnew+B*u(index); k4=step*F; x=x+(1/6)*(k1+2*k2+2*k3+k4); % Calculated next state t=t+step; % Increments time index=index+1; % Increments the indexing term so that % index=1 corresponds to time t=0. end % This end statement goes with the first "while" statement % in the program % % Next, we provide plots of the input and output of the plant % and the output of the reference model % subplot(211) plot(time,y,'-',time,ym,'--') xlabel('Time (sec)') title('Output of plant (solid) and reference model (dashed)') subplot(212) plot(time,u) xlabel('Time (sec)') title('Output of fuzzy controller (input to the plant)') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % End of program % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %***************************************************** >>--->>CUT-HERE>>--------------------------------- %***************************************************** /* E6.1 */ /*********************************************************************** * C version: by Min-Hsiung Hung 6/16/97 * Reprogrammed from the Matlab version: * Brian Klinehoffer and Kevin Passino, 5/12/97. * ************************************************************************ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Fuzzy Model Reference Learning Control System Simulator %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Matlab version: Brian Klinehoffer and Kevin Passino % Version: 5/12/97 % % This program is used to simulate the "fuzzy model reference % learning controller" (FMRLC) for a simple plant. The objective % here is to illustrate the basic issues in how to simulate this % adaptive fuzzy controller and to show how it will perform under % various conditions. % % The plant is a second order system in the form of % k_p/(s^2 + 2*zeta_p*w_p*s + w_p^2) where k_p is the gain in the % numerator, zeta_p is the damping ratio and w_p is the undamped % natural frequency. We encourage the user to change these parameters % (e.g., even on-line) to see how the FMRLC can adapt to the % changes in the plant that such a variation represents. We use % the simple linear system for the sake of illustration. We have % found the FMRLC particularly well-suited for a variety of nonlinear % plants and it is for these that the FMRLC shows some distinct % advantages over conventional methods (for some applications). % Suppose that we name the input to the plant u and its output y. % % The goal of the FMRLC is to make the closed-loop system behave like % a given "reference model." Here, our reference model is a % first order system that represents how we would like the system % to behave. For instance, the plant may have unknown pole positions % with a low zeta_p value so that the plant is oscillatory. The % reference model could specify that we would like the pole of % the closed-loop system to be further out into the left-half plane % and have a higher zeta so that it will be faster. The input to % the reference model is the reference input r and its output is ym. % % There are four components to the FMRLC: the plant, reference model, % fuzzy controller that will be tuned, and "fuzzy inverse % model." We will not fully explain the rationale for the choice of % the parameters for each of these in this program. Our objective % here is simply to show how to program the FMRLC once you have a % basic understanding of its operation. % % The fuzzy controller here has two inputs and one output. One input % is the error e=r-y and the other is the change in error % c=(e(kT)-e(kt-T))/T where T=step (in our program). The fuzzy % controller uses 11 uniformly distributed triangular membership % functions for each of its input universes of discourse. We use % "linguistic-numeric" indices of 1,2,3,...,11 to name the membership % functions (similarly for the fuzzy inverse model). % The fuzzy controller uses minimum to represent the premise % and implication, and COG defuzzification. The output membership % function centers of the fuzzy controller are what is tuned by the % learning mechanism of the FMRLC. Clearly, all these choices for % the structure and parameters are simply for illustrative purposes; % other choices may be more appropriate for other plants and control % objectives. % % The learning mechanism consists of two parts: the fuzzy inverse % model and the knowledge-base modifier. The fuzzy inverse model % is simply a fuzzy controller that operates in the adaptation loop. % Its two inputs are ye=ym - y and yc=(ye(kT)-ye(kT-T))/T (again % T=step). It uses 11 uniformly distributed triangular membership % functions for each of its input universes of discourse. It uses % triangular shaped output membership functions. It uses minimum % to represent the premise and implication, and COG defuzzification. % The output of the fuzzy inverse model is p and it indicates how % much to move the centers of the output membership functions. % The knowledge-base modifier is the mechanism that actually moves % the fuzzy controller output membership function centers. Here, % we use the fact that either 1, 2, or 4 rules of the fuzzy % controller will be on at any one time. Our update formula is % c_i(kT)=c_i(kT-dT)+p(kT) where c_i is the output center for % the ith rule and the index i is selected from the set of rules % that were on d steps before (at most 4 rules) (we use different % notation in the program below). In this way the rules that % contributed to the current error y_e are the ones that % are modified by the learning mechanism. It is important to note % that while for the simple linear system we do not design a very % innovative fuzzy inverse model, real-world applications will % dictate the need for much more careful tuning/design of the % rule-base and membership functions for the inverse model. The % one we have chosen works adequately for this simple problem. % Next, note that in the code we have taken certain shortcuts % soley for illustrative purposes. For instance, the code for the % fuzzy controller that is tuned is very similar to that of the % fuzzy inverse model; while this works for this simple application % clearly it will not always be possible to use such similar % fuzzy systems (e.g., ones that have the same number of inputs, % rules, and membership functions). This completes the description % of the FMRLC. % % The FMRLC is simulated with a 4th order Runge-Kutta technique. One % must be careful in the choice of the integration step size. For % instance, if you change the a_p parameter to be very large, which % corresponds to a fast system, then you will need to make the value % of the integration step size smaller so that the simulation is % accurate. This completes our brief description of how this program % operates; however, comments are included throughout the program to % further clarify its operation. If you want more details on how the % FMRLC works see the following references: % % Layne J.R., Passino K.M., "Fuzzy Model Reference Learning Control % for Cargo Ship Steering", IEEE Control Systems Magazine, Vol. 13, % No. 6, pp. 23-34, Dec. 1993. % % Layne J.R., Passino K.M., "Fuzzy Model Reference Learning Control," % To appear in the Journal of Intelligent and Fuzzy Systems, 1996. % % Layne J.R., Passino K.M., Yurkovich S., "Fuzzy Learning Control % for Anti-Skid Braking Systems", IEEE Trans. on Control Systems % Technology, Vol. 1, No. 2, pp. 122-129, June 1993. % % Kwong W.A., Passino K.M., Lauknonen E.G., Yurkovich S., ³Expert % Supervision of Fuzzy Learning Systems for Fault Tolerant Aircraft % Control², Special Issue on Fuzzy Logic in Engineering Applications, % Proceedings of the Inst. of Electrical and Electronics Engineers, % Vol. 83, No.3, pp. 466-483, March 1995. % % Moudgal V.G., Kwong W.A., Passino K.M., and Yurkovich S., ³Fuzzy % Learning Control for a Flexible-Link Robot², IEEE Transactions on % Fuzzy Systems, Vol. 3, No. 2, pp. 199-210, May 1995. % % Kwong W.A., Passino K.M., "Dynamically Focused Fuzzy Learning % Control", IEEE Trans. on Systems, Man, and Cybernetics, Vol. 26, % No. 1, pp. 53-74, Feb. 1996. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ #include #include #include #define MAX(A,B) ((A) > (B) ? (A) : (B)) #define MIN(A,B) ((A) < (B) ? (A) : (B)) /* function protocals */ void ec_centers(double *ce, double ge, double *cc, double gc); // setting centers of mfs for e & c void fc_rules(double fuzzyrules[][11], double gu, double gf); // establish the rule table for the fuzzy controller void yeyc_centers(double *cye, double gye, double *cyc, double gyc); // setting centers of mfs for ye & yc void finvm_rules(double inverrules[][11], double gp); // establish the rule table for the fuzzy inverse model void main(void) { //******************* variable declarations ***************************************** // Initialize variables double eold=0.0, // Intial condition used to calculate c yeold=0.0, // Intial condition used to calculate yc ymold=0.0; // Initial condition for the first order reference model // Next, initialize parameters for the fuzzy controller int nume=11, // Number of input membership functions for the e // universe of discourse numc=11; // Number of input membership functions for the c // universe of discourse double ge=1.0/2.0, gc=1.0/2.0, gu=5.0; // Scaling gains for tuning membership functions for // universes of discourse for e, c and u respectively // These are tuned to improve the performance of the FMRLC double we=0.2*(1.0/ge); // we is half the width of the triangular input membership // function bases (note that if you change ge, the base width // will correspondingly change so that we always end // up with uniformly distributed input membership functions) // Note that if you change nume you will need to adjust the // "0.2" factor if you want membership functions that // overlap in the same way. double wc=0.2*(1.0/gc); // Similar to we but for the c universe of discourse double base=0.4*gu; // Base width of output membership fuctions of the fuzzy // controller // Place centers of membership functions of the fuzzy controller: // Centers of input membership functions for the e universe of // discourse for of fuzzy controller (a vector of centers) double ce[11], //The values will be set in the "ec_centers()" routine // Centers of input membership functions for the c universe of // discourse for of fuzzy controller (a vector of centers) cc[11]; //The values will be set in the "ec_centers()" routine // This next matrix determines the rules of the fuzzy controller. // The entries are the centers of the output membership functions. // This choice represents our "best guess" on how to synthesize // the fuzzy controller for the nominal system (with a little // work you could certainly do better). Notice the regularity // of the pattern of rules (basiscally we are using a type of // saturated index adding). Notice that it is scaled by gu, the // output scaling factor, since it is a normalized rule-base // // Also, there is the parameter gf that is either =0 or =1. It // mulitplies the rule-base and is =1 if you want the FMRLC to be // initialized with a decent guess at what the fuzzy controller // that it is going to tune should be. If gf=0 then the rule-base // has 121 rules with output membership functions that are triangular // with base widths base=0.4*gu and centers at zero. Hence, // if gf=0 this represents that the fuzzy controller does not // know very much at all about how to control the plant initially. // As the FMRLC learns it will fill in the rule-base to // try to achieve good behavior. double gf=0.0; double fuzzyrules[11][11]; //The values will be set in the "fc_rules()" routine // Next, we define some parameters for the fuzzy inverse model double gye=1.0/2.0, gyc=1.0/2.0; // Scaling gains for the error and change in error for // the inverse model // These are tuned to improve the performance of the FMRLC double gp=0.2; // Scaling gain for the output of inverse model. If you let gp=0 then // you turn off the learning mechanism and hence the FMRLC reduces to // a direct fuzzy controller that is not tuned. Hence, from this // program you can also see how to code a non-adaptive fuzzy controller. // A reasonable value for it is gp=.2 but you should try other values. int numye=11, // Number of input membership functions for the ye // universe of discourse numyc=11; // Number of input membership functions for the yc // universe of discourse double wye=0.2*(1.0/gye), // Sets the width of the membership functions for // ye from center to extremes wyc=0.2*(1.0/gyc), // Sets the width of the membership functions for // yc from center to extremes invbase=0.4*gp; // Sets the base of the output membership functions // for the inverse model // Place centers of inverse model membership functions // For error input for learning Mechanism double cye[11], //The values will be set in the "yeyc_centers()" routine // For change in error input for learning mechanism cyc[11]; //The values will be set in the "yeyc_centers()" routine // The next matrix contains the rule-base matrix for the fuzzy // inverse model. Notice that for simplicity we choose it to have // the same structure as the rule-base for the fuzzy controller. // While this will work for the control of the simple first order // linear system for many nonlinear systems a different structure // will be needed for the rule-base. Again, the entries are // the centers of the output membership functions, but now for // the fuzzy inverse model. double inverrules[11][11]; //The values will be set in the "finvm_rules()" routine // Next, we set up some parameters/variables for the // knowledge-base modifier int d=1; // This sets the number of steps the knowledge-base modifier looks // back in time. For this program it must be an integer // less than or equal to 10 (but this is easy to make larger) // The next four vectors are used to store the information about // which rules were on 1 step in the past, 2 steps in the past, ...., // 10 steps in the past (so that picking 0<= d <= 10 can be used). int meme_int[10]={0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; // sets up the vector to store up to 10 values of e_int int meme_count[10]={0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; // sets up the vector to store up to 10 values of e_count int memc_int[10]={0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; // sets up the vector to store up to 10 values of c_int int memc_count[10]={0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; // sets up the vector to store up to 10 values of c_count // // Next, we intialize the simulation of the closed-loop system. // #define k_p 1.0 // The numerator of the plant. Change this value to study // the ability of the FMRLC to control other plants. Also, // you can make this a time-varying parameter. #define zeta_p 0.707 // Damping ratio for the second order plant (could change this // to see how the system will adapt to it) #define w_p 1.0 // Undamped natural frequency for the plant (could change this // to see how the system will adapt to it) #define k_r 1.0 // The numerator of the reference model. Change this value to study // the ability of the FMRLC to meet other performance specifications. #define a_r 1.0 // The value of -a_r is the pole position for the reference model. // Change this value to study the ability of the FMRLC to meet // other performance specifications (e.g., a faster response). double t=0.0; // Reset time to zero int index=0; // This is time's index (not time, its index). double tstop=64.0; // Stopping time for the simulation (in seconds) double step=0.01; // Integration step size double x[2]={0,0}; // Intial condition on state of the plant // Need a state space representation for the plant. Since our // plant is linear we use the standard form of xdot=Ax+Bu, y=Cx+Du // Matrix A of state space representation of plant double A[2][2]={ {0.0 , 1.0 }, {-1*w_p*w_p, -2*zeta_p*w_p} }; // Matrix A of state space representation of plant double B[2]={0.0, 1.0}; // Matrix B of state space representation of plant double C[2]={k_p, 0.0}; // Matrix C of state space representation of plant // Declarations of other variables used in the program int i,j,k, l; // used in loop double y, r, ym, u; // recorded variables: y : output of the plant // r : reference input // ym: output of the reference model // u : plant input double rold=0.0, e, c, num, den, prem, ye, yc, invnum, invden, p; int e_count , c_count, ye_count, yc_count, e_int, c_int, ye_int, yc_int; double mfe[11], mfc[11], mfye[11], mfyc[11]; double F[2], k1[2], k2[2], k3[2],k4[2], xnew[2]; // used in Runge-Kutta equations FILE *myfile; //***************** end of variable declarations ************************************* // Excuting the initializing routines: ec_centers(ce,ge,cc,gc); // setting centers of mfs for e & c fc_rules(fuzzyrules, gu, gf); // establish the rule table for the fuzzy controller yeyc_centers(cye,gye,cyc,gyc); // setting centers of mfs for ye & yc finvm_rules(inverrules, gp); // establish the rule table for the fuzzy inverse model // check the initialization //----- ce & ge ---------------- printf("ce=[ "); for (i=0; i<=10; i++) printf("%f ", ce[i]); printf(" ], ge=%f\n\n",ge); //----- cc & gc ---------------- printf("cc=[ "); for (i=0; i<=10; i++) printf("%f ",cc[i]); printf(" ], gc=%f\n\n", gc); //----- cye & gye ---------------- printf("cye=[ "); for (i=0; i<=10; i++) printf("%f ", cye[i]); printf(" ], gye=%f\n\n",gye); //----- cyc & gyc ---------------- printf("cyc=[ "); for (i=0; i<=10; i++) printf("%f ",cyc[i]); printf(" ], gyc=%f\n\n", gyc); //------- gu, gf, gp ---------------- printf("gu=%f, gf=%f, gp=%f\n\n", gu, gf, gp); //------- A, B, C -------------------- for (i=0; i<=1; i++) for (j=0; j<=1; j++) printf("A[%d][%d]=%f\n", i, j, A[i][j]); printf("\n"); printf("B[0]=%f, B[1]=%f\n\n", B[0], B[1]); printf("C[0]=%f, C[1]=%f\n\n", C[0], C[1]); //--------- fuzzyrules -------------------------- printf("fuzzyrules=\n"); for (i=0; i<=10; i++) { for (j=0; j<=9; j++) printf("%f ", fuzzyrules[i][j]); printf("%f\n",fuzzyrules[i][10]); } //--------- inverrules -------------------------- printf("inverrules=\n"); for (i=0; i<=10; i++) { for (j=0; j<=9; j++) printf("%f ", inverrules[i][j]); printf("%f\n", inverrules[i][10]); } // // Next, we start the simulation of the system. This is the main // loop for the simulation of the FMRLC. // myfile=fopen("data.dat","w"); // open a file to store data while (t <= tstop ) { y=C[0]*x[0]+C[1]*x[1]; // Output of the plant // Next, we define the reference input r as a sine wave r=sin(0.6*t); // The following sets up the reference model. Just for the sake of // illustration we use a discrete-time reference model. If you use // a continuous-time reference model you will have to augment // the state of the closed-loop system with the state(s) of the // reference model and hence update the state in the Runge-Kutta // equations above. // // For the reference model we use a first order transfer function // k_r/(s+a_r) but we use the bilinear transformation where we // replace s by (2/step)(z-1)/(z+1), then find the z-domain // representation of the reference model, then convert this // to a difference equation: ym=(1.0/(2.0+a_r*step))*((2.0-a_r*step)*ymold+k_r*step*(r+rold)); ymold=ym; rold=r; // This saves the past value of the ym and r so that we can use it // the next time around the loop // Now that we have simulated the next step for the plant and reference // model we will focus on the two fuzzy components. // First, for the given fuzzy controller inputs we determine // the extent at which the error membership functions // of the fuzzy controller are on (this is the fuzzification part). c_count=0; e_count=0; // These are used to count the number of // non-zero mf certainities of e and c e=r-y; // Calculates the error input for the fuzzy controller c=(e-eold)/step; // Calculates the change in error input for the fuzzy controller eold=e; // Saves the past value of e for use in the next time through the // loop // The following if-then structure fills the vector mfe // with the certainty of each membership fucntion of e for the // current input e if (e<=ce[0]) // Takes care of saturation of the left-most { // membership function mfe[0]=1.0; // i.e., the only one on is the left-most one for (i=1; i<=10; i++) mfe[i]=0.0; e_count=e_count+1; e_int=0; } else { if (e>=ce[nume-1]) // Takes care of saturation of the right-most mf { for (i=0; i<=9; i++) mfe[i]=0.0; mfe[10]=1.0; // One mf on, it is the right-most one e_count=e_count+1; e_int=nume-1; } else // In this case the input is on the middle part of the { // universe of discourse for e // Next, we are going to cycle through the mfs to // find all that are on for (i=0; i<=nume-1; i++) { if (e<=ce[i]) { mfe[i]=MAX(0, (1.0+(e-ce[i])/we) ); // In this case the input is to the // left of the center ce(i) and we compute // the value of the mf centered at ce(i) // for this input e if (mfe[i]!=0) { // If the certainty is not equal to zero then say // that have one mf on by incrementing our count e_count=e_count+1; e_int=i; // This term holds the index last entry // with a non-zero term } } else { mfe[i]=MAX(0,(1.0+(ce[i]-e)/we) ); // In this case the input is to the // right of the center ce(i) if (mfe[i]!=0) { e_count=e_count+1; e_int=i; // This term holds the index of the // last entry with a non-zero term } } // end of else } // end of for } // end of else } // end of else // Next we will save the number of mfs that are on and the pointer // e_int as to which rules were on. This vector of length // 10 saves the last 10 values of e_count and e_int as time // progresses (hence, it is a moving window). // These will be used by the FMRLC knowledge-base modifier. for (i=0; i<=8; i++) { meme_count[i+1]=meme_count[i]; meme_int[i+1]=meme_int[i]; } meme_count[0]=e_count; meme_int[0]=e_int; // The following if-then structure fills the vector mfc with the // certainty of each membership fucntion of c for the current // value of c (to understand this part of the code see the above // similar code for computing mfe) if (c<=cc[0]) // Takes care of saturation of the left-most { // membership function mfc[0]=1.0; // i.e., the only one on is the left-most one for (i=1; i<=10; i++) mfc[i]=0.0; c_count=c_count+1; c_int=0; } else { if (c>=cc[numc-1]) // Takes care of saturation of the right-most mf { for (i=0; i<=9; i++) mfc[i]=0.0; mfc[10]=1.0; // One mf on, it is the right-most one c_count=c_count+1; c_int=numc-1; } else // In this case the input is on the middle part of the { // universe of discourse for c // Next, we are going to cycle through the mfs to // find all that are on for (i=0; i<=numc-1; i++) { if (c<=cc[i]) { mfc[i]=MAX(0, (1.0+(c-cc[i])/wc) ); // In this case the input is to the // left of the center cc(i) and we compute // the value of the mf centered at cc(i) // for this input c if (mfc[i]!=0) { // If the certainty is not equal to zero then say // that have one mf on by incrementing our count c_count=c_count+1; c_int=i; // This term holds the index last entry // with a non-zero term } } else { mfc[i]=MAX(0,(1.0+(cc[i]-c)/wc) ); // In this case the input is to the // right of the center cc(i) if (mfc[i]!=0) { c_count=c_count+1; c_int=i; // This term holds the index of the // last entry with a non-zero term } } // end of else } // end of for } // end of else } // end of else // Next we will save the number of mfs that are on and the pointer // c_int as to which rules were on. This vector of length 10 // saves the last 10 values of c_count and e_int as time progresses // (hence, it is a moving window). These will be used by the FMRLC // knowledge-base modifier. for (i=0; i<=8; i++) { memc_count[i+1]=memc_count[i]; memc_int[i+1]=memc_int[i]; } memc_count[0]=c_count; memc_int[0]=c_int; // These for loops calculate the crisp output using only the non- // zero premise of error,e, and change in error,c. This cuts down // computation time since we will only compute // the contribution from the rules that are on (i.e., a maximum of // four rules for our case). Minimum is used for the premise // and implication. num=0.0; den=0.0; for (k=(e_int-e_count+1); k<=e_int; k++) { // Scan over e indices of mfs that are on for (l=(c_int-c_count+1); l<=c_int; l++) { // Scan over c indices of mfs that are on prem=MIN(mfe[k], mfc[l]); // Value of premise membership function // This next calculation of num adds up the numerator for the // defuzzification formula. fuzzyrules(k,l) is the output center // for the rule.base*(prem-(prem)^2/2 is the area of a symmetric // triangle with base width "base" and that is chopped off at // a height of prem (since we use minimum to represent the // implication). Computation of den is similar but without // fuzzyrules(k,l). num=num+fuzzyrules[k][l]*base*(prem-(prem*prem)/2.0); den=den+base*(prem-(prem*prem)/2.0); } } u=num/den; // Crisp output of fuzzy controller that is the input // to the plant // Next, we perform computations for the fuzzy inverse model. ye=ym-y; // Calculates ye yc=(ye-yeold)/step; // Calculates yc yeold=ye; // Saves the value of ye for use the next time ye_count=0; yc_count=0; // Counts the non-zero certainties // of ye and yc similar to how we did // for the fuzzy controller // The following if-then structure fills the vector mfye with the // certainty of each membership fucntion of ye (similar to the // fuzzy controller). Notice that we use the same number of // input membership functions as in the fuzzy controller - // just for convenience - it does not have to be this way if (ye<=cye[0]) // Takes care of saturation of the left-most { // membership function mfye[0]=1.0; // i.e., the only one on is the left-most one for (i=1; i<=10; i++) mfye[i]=0.0; ye_count=ye_count+1; ye_int=0; } else { if (ye>=cye[numye-1]) // Takes care of saturation of the right-most mf { for (i=0; i<=9; i++) mfye[i]=0.0; mfye[10]=1.0; // One mf on, it is the right-most one ye_count=ye_count+1; ye_int=numye-1; } else // In this case the input is on the middle part of the { // universe of discourse for ye // Next, we are going to cycle through the mfs to // find all that are on for (i=0; i<=numye-1; i++) { if (ye<=cye[i]) { mfye[i]=MAX(0, (1.0+(ye-cye[i])/wye) ); // In this case the input is to the // left of the center cye(i) and we compute // the value of the mf centered at cye(i) // for this input ye if (mfye[i]!=0) { // If the certainty is not equal to zero then say // that have one mf on by incrementing our count ye_count=ye_count+1; ye_int=i; // This term holds the index last entry // with a non-zero term } } else { mfye[i]=MAX(0,(1.0+(cye[i]-ye)/wye) ); // In this case the input is to the // right of the center cye(i) if (mfye[i]!=0) { ye_count=ye_count+1; ye_int=i; // This term holds the index of the // last entry with a non-zero term } } // end of else } // end of for } // end of else } // end of else // The following if-then structure fills the vector mfyc with the // certainty of each membership fucntion of yc if (yc<=cyc[0]) // Takes care of saturation of the left-most { // membership function mfyc[0]=1.0; // i.e., the only one on is the left-most one for (i=1; i<=10; i++) mfyc[i]=0.0; yc_count=yc_count+1; yc_int=0; } else { if (yc>=cyc[numyc-1]) // Takes care of saturation of the right-most mf { for (i=0; i<=9; i++) mfyc[i]=0.0; mfyc[10]=1.0; // One mf on, it is the right-most one yc_count=yc_count+1; yc_int=numyc-1; } else // In this case the input is on the middle part of the { // universe of discourse for yc // Next, we are going to cycle through the mfs to // find all that are on for (i=0; i<=numyc-1; i++) { if (yc<=cyc[i]) { mfyc[i]=MAX(0, (1.0+(yc-cyc[i])/wyc) ); // In this case the input is to the // left of the center cyc(i) and we compute // the value of the mf centered at cyc(i) // for this input yc if (mfyc[i]!=0) { // If the certainty is not equal to zero then say // that have one mf on by incrementing our count yc_count=yc_count+1; yc_int=i; // This term holds the index last entry // with a non-zero term } } else { mfyc[i]=MAX(0,(1.0+(cyc[i]-yc)/wyc) ); // In this case the input is to the // right of the center cyc(i) if (mfyc[i]!=0) { yc_count=yc_count+1; yc_int=i; // This term holds the index of the // last entry with a non-zero term } } // end of else } // end of for } // end of else } // end of else // These for loops calculate the crisp output of the inverse model // using only the non-zero premise of error,ye, and change in // error,yc. This cuts down computation time (similar to the // fuzzy controller). Minimum is used for the premise and the // implication. To understand this part of the code see the // similar code for the fuzzy controller. invnum=0.0; invden=0.0; for (k=(ye_int-ye_count+1); k<=ye_int; k++) { for (l=(yc_int-yc_count+1); l<=yc_int; l++) { prem=MIN(mfye[k], mfyc[l]); invnum=invnum+inverrules[k][l]*invbase*(prem-(prem*prem)/2.0); invden=invden+invbase*(prem-(prem*prem)/2.0); } } // Next we compute the output of the fuzzy inverse model. // If you want to let gp=0 to test the fuzzy controller by itself // then you will have invden=0 and you will not be able to compute p. // To make it possible to let gp=0 we put in the following // if-then rule. if (gp==0) p=0; else p=invnum/invden; // Crisp output of inverse model //if (abs(p)<.05) p=0; // robustification term // Next we modify the centers of ouput membership functions // associated with the rules that were on d steps ago // The indexing sheme is similar to the one used in the // computation of the outputs of the fuzzy controller and // fuzzy inverse model. for (k=(meme_int[d]-meme_count[d]+1); k<=meme_int[d]; k++) for (l=(memc_int[d]-memc_count[d]+1); l<=memc_int[d]; l++) fuzzyrules[k][l]=fuzzyrules[k][l]+p; // Next, the Runge-Kutta equations to find the next state. // k1 F[0]=A[0][0]*x[0]+A[0][1]*x[1]+B[0]*u; F[1]=A[1][0]*x[0]+A[1][1]*x[1]+B[1]*u; k1[0]=step*F[0]; k1[1]=step*F[1]; xnew[0]=x[0]+k1[0]/2.0; xnew[1]=x[1]+k1[1]/2.0; // k2 F[0]=A[0][0]*xnew[0]+A[0][1]*xnew[1]+B[0]*u; F[1]=A[1][0]*xnew[0]+A[1][1]*xnew[1]+B[1]*u; k2[0]=step*F[0]; k2[1]=step*F[1]; xnew[0]=x[0]+k2[0]/2.0; xnew[1]=x[1]+k2[1]/2.0; // k3 F[0]=A[0][0]*xnew[0]+A[0][1]*xnew[1]+B[0]*u; F[1]=A[1][0]*xnew[0]+A[1][1]*xnew[1]+B[1]*u; k3[0]=step*F[0]; k3[1]=step*F[1]; xnew[0]=x[0]+k3[0]; xnew[1]=x[1]+k3[1]; // k4 F[0]=A[0][0]*xnew[0]+A[0][1]*xnew[1]+B[0]*u; F[1]=A[1][0]*xnew[0]+A[1][1]*xnew[1]+B[1]*u; k4[0]=step*F[0]; k4[1]=step*F[1]; // update x x[0]=x[0]+(1.0/6.0)*(k1[0]+2.0*k2[0]+2.0*k3[0]+k4[0]); // Calculated next state x[1]=x[1]+(1.0/6.0)*(k1[1]+2.0*k2[1]+2.0*k3[1]+k4[1]); // Calculated next state // save data fprintf(myfile, "%f %f %f %f %f\n", t, y, ym, u, r); t=t+step; // Increments time index=index+1; // Increments the indexing term so that // index=0 corresponds to time t=0. } // This end statement goes with the first "while" statement in the program fclose(myfile); // close the data file // save the rule-base matrix of the fuzzy controller into a disk file myfile=fopen("fuzrules.dat","w"); for (i=0; i<=10; i++) { for (j=0; j<=9; j++) fprintf(myfile,"%f ", fuzzyrules[i][j]); fprintf(myfile,"%f\n",fuzzyrules[i][10]); } fclose(myfile); // save the rule-base matrix of the fuzzy inverse model into a disk file myfile=fopen("invrules.dat","w"); for (i=0; i<=10; i++) { for (j=0; j<=10; j++) fprintf(myfile,"%f ", inverrules[i][j]); fprintf(myfile,"%f\n",inverrules[i][10]); } fclose(myfile); } // This end the main function //**************** The followings are subroutines ********************************* // Setting the center values of the input membership functions for the e & c // universe of discourse void ec_centers(double *ce, double ge, double *cc, double gc) { int i; // Centers of input membership functions for the e universe of // discourse for of fuzzy controller (a vector of centers) for (i=0; i<=10; i++) ce[i]=(-1.0+0.2*i)*(1.0/ge); // Centers of input membership functions for the c universe of // discourse for of fuzzy controller (a vector of centers) for (i=0; i<=10; i++) cc[i]=(-1.0+0.2*i)*(1.0/gc); } // Setting the initial values of the rule-base matrix for the fuzzy controller. void fc_rules(double fuzzyrules[][11], double gu, double gf) { // This next matrix determines the rules of the fuzzy controller. // The entries are the centers of the output membership functions. // This choice represents our "best guess" on how to synthesize // the fuzzy controller for the nominal system (with a little // work you could certainly do better). Notice the regularity // of the pattern of rules (basiscally we are using a type of // saturated index adding). Notice that it is scaled by gu, the // output scaling factor, since it is a normalized rule-base // // Also, there is the parameter gf that is either =0 or =1. It // mulitplies the rule-base and is =1 if you want the FMRLC to be // initialized with a decent guess at what the fuzzy controller // that it is going to tune should be. If gf=0 then the rule-base // has 121 rules with output membership functions that are triangular // with base widths base=0.4*gu and centers at zero. Hence, // if gf=0 this represents that the fuzzy controller does not // know very much at all about how to control the plant initially. // As the FMRLC learns it will fill in the rule-base to // try to achieve good behavior. int i,j; double temp[11][11]= { {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0}, {-1.0, -1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2}, {-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4}, {-1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6}, {-1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8}, {-1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0}, {-0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0}, {-0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0}, {-0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0}, {-0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0, 1.0}, { 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0} }; for (i=0;i<=10;i++) for (j=0; j<=10; j++) fuzzyrules[i][j]=temp[i][j]*gu*gf; } // Setting the center values of the input membership functions for the ye & yc // universe of discourse void yeyc_centers(double *cye, double gye, double *cyc, double gyc) { int i; // Centers of input membership functions for the e universe of // discourse for of fuzzy controller (a vector of centers) for (i=0; i<=10; i++) cye[i]=(-1.0+0.2*i)*(1.0/gye); // Centers of input membership functions for the c universe of // discourse for of fuzzy controller (a vector of centers) for (i=0; i<=10; i++) cyc[i]=(-1.0+0.2*i)*(1.0/gyc); } // Setting the values of the rule-base matrix for the fuzzy inverse model. void finvm_rules(double inverrules[][11], double gp) { // The next matrix contains the rule-base matrix for the fuzzy // inverse model. Notice that for simplicity we choose it to have // the same structure as the rule-base for the fuzzy controller. // While this will work for the control of the simple first order // linear system for many nonlinear systems a different structure // will be needed for the rule-base. Again, the entries are // the centers of the output membership functions, but now for // the fuzzy inverse model. int i,j; double temp[11][11]= { {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0}, {-1.0, -1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2}, {-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4}, {-1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6}, {-1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8}, {-1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0}, {-0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0}, {-0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0}, {-0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0}, {-0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0, 1.0}, { 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0} }; for (i=0;i<=10;i++) for (j=0; j<=10; j++) inverrules[i][j]=temp[i][j]*gp; } /* %************ Drawing Routine (MATLAB) for the FMRLC simulation *********** clear; clc; load data.dat % load data: t, y, ym, u, r load fuzrules.dat % load rule_base matrix of the fuzzy controller load invrules.dat % load rule_base matrix of the fuzzy inverse model t=data(:,1); y=data(:,2); ym=data(:,3); u=data(:,4); r=data(:,5); subplot(211) plot(t,y,'-',t,ym,'--') xlabel('Time (sec)') title('Output of plant (solid) and reference model (dashed)') subplot(212) plot(t,u) xlabel('Time (sec)') title('Output of fuzzy controller (input to the plant)') */ %***************************************************** >>--->>CUT-HERE>>--------------------------------- %***************************************************** % END OF FILE