diff options
Diffstat (limited to 'Process_Dynamics_And_Controls_by_D_E_Seborg')
16 files changed, 5040 insertions, 0 deletions
diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/10-Dynamic_behavior.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/10-Dynamic_behavior.ipynb new file mode 100644 index 0000000..99f82f7 --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/10-Dynamic_behavior.ipynb @@ -0,0 +1,348 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 10: Dynamic behavior" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 10.10: Routh_Array_1.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 10.10\n", +"disp('Example 10.10')\n", +"s=%s;\n", +"Gp=1/(5*s+1);\n", +"Gm=1/(s+1);\n", +"Gv=1/(2*s+1);\n", +"Ys=Gv*Gp*Gm\n", +"Routh=routh_t(Ys,poly(0,'Kc')); // produces routh table for polynomial 1+Kc*Ys\n", +"disp(Routh)\n", +"K1=roots(numer(Routh(3,1)));\n", +"K2=roots(numer(Routh(4,1)));\n", +"mprintf('K lies between %f and %f for system to be stable', K2,K1)" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 10.11: Routh_Array_2.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 10.11\n", +"disp('Example 10.11')\n", +"Kc=poly(0,'Kc');//defining a polynomial variable\n", +"a2=2.5;a1=5.5-Kc;a0=1+2*Kc;//a# are coefficients\n", +"b1=(a1*a0-a2*0)/a1;\n", +"mprintf('Routh Array is')\n", +"A=[a2 a0;a1 0;b1 0]\n", +"disp(A)\n", +"mprintf('\n All entries in first column should be positive')\n", +"Kc_up=roots(a1);//upper limit for Kc by solving third row column 1 of array\n", +"b1=numer(b1);//This is done to extract the numerator from rational c1\n", +"//without extracting numerator we cannot find roots using 'roots' function\n", +"Kc_ul=roots(b1);//lower limit for Kc\n", +"mprintf('\n\nThe values of Kc for system to be stable are \n %f<Kc<%f',Kc_ul,Kc_up);" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 10.12: Direct_substitution_to_find_stability.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 10.12\n", +"disp('Example 10.12')\n", +"w=poly(0,'w')\n", +"s=%i*w; //j or iota is i\n", +"G=10*s^3+17*s^2+8*s+1;//Kc has been removed here because in a single expression\n", +"//two polynomials are not allowed\n", +"w=roots(imag(G));\n", +"p=-real(G)//Real part of G\n", +"Kc=horner(p,w)\n", +"mprintf('The values outside which system is unstable \nare %f and %f',Kc(1),Kc(3))" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 10.13: Root_Locus.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 10.13\n", +"disp('Example 10.13')\n", +"s=%s;\n", +"G=4/((s+1)*(s+2)*(s+3));\n", +"G=syslin('c',G);\n", +"[ki,s_i]=kpure(G);\n", +"evans(G,ki*1.5); // plots for until K = 1.5*ki\n", +"//disp(G,'For G=');disp(ki,'K=')\n", +"disp(ki,'Max value of k for which we have closed loop stability',G,'For G=')\n", +"xtitle('$G(s)=\frac{4}{(s+1)(s+2)(s+3)}$')\n", +"sgrid();" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 10.14: Transient_response_from_root_locus.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 10.14\n", +"disp('Example 10.14')\n", +"s=%s;\n", +"G=4/((s+1)*(s+2)*(s+3));\n", +"K=10; //given in question\n", +"p=1+K*G;//characteristic equation\n", +"q=roots(numer(p));\n", +"q_abs=abs(q);\n", +"q_real=real(q);\n", +"q_imag=imag(q);\n", +"d=q_abs(2);\n", +"psi=%pi-acos(q_real./q_abs);//angle in radians\n", +"tau=1/d;\n", +"eta=cos(psi)\n", +"mprintf('\nd=%f\npsi=%f degrees\ntau=%f time units\neta=%f',d,psi(2)*180/%pi,tau,eta(2))\n", +"mprintf('\n\nPlease note that the answers given in book are wrong')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 10.2: Set_point_response_of_level_control_system.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 10.2\n", +"disp('Example 10.2')\n", +"A=%pi*0.5^2; //Square meters\n", +"R=6.37;\n", +"Kp=R//min/sq.m=R\n", +"tau=R*A;\n", +"Km=100/2; //% per meter\n", +"l=0.5;\n", +"q=0.2*30^(l-1);\n", +"dq_dl=0.2*log(30)*30^(l-1); //cu.meter/min Eqn 10-48\n", +"Kip=(15-3)/100;//psi/%\n", +"dl_dpt=(1-0)/(15-3); //psi^-1\n", +"Kv=dq_dl*dl_dpt //Eqn 10-50\n", +"Kc=[4 8 20]'; //different values of Kc that we have to try\n", +"K_OL=Kc*Kv*Kp*Km*Kip;//Open loop gain Eqn 10-40\n", +"K1=K_OL./(1+K_OL);//Eqn 10-38\n", +"tau1=tau./(1+K_OL); //Eqn 10-39\n", +"//Now we simulate the close loop process for these different values of K1 and tau1\n", +"s=%s;\n", +"G=K1./(tau1*s+1);\n", +"G=syslin('c',G);\n", +"t=[0:0.1:10]'; //time in minutes\n", +"hdash=csim('step',t,G)'; \n", +"plot2d(t,hdash,rect=[0,0,10,1.25])\n", +"xgrid()\n", +"xtitle('Ex-10.2','Time(min)','h(t)');\n", +"a=legend('Kc=4','Kc=8','Kc=20',position=4);\n", +"a.font_size=2;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 10.3: Disturbance_response_of_level_control_system.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 10.3\n", +"disp('Example 10.3')\n", +"//This example draws upon the calculations of Ex 10.2 and hence it has been \n", +"//reproduced\n", +"A=%pi*0.5^2; //Square meters\n", +"R=6.37;\n", +"Kp=R//min/sq.m=R\n", +"tau=R*A;\n", +"Km=100/2 //% per meter\n", +"l=0.5;\n", +"q=0.2*30^(l-1);\n", +"dq_dl=0.2*log(30)*30^(l-1); //cu.meter/min Eqn 10-48\n", +"Kip=(15-3)/100;//psi/%\n", +"dl_dpt=(1-0)/(15-3); //psi^-1\n", +"Kv=dq_dl*dl_dpt //Eqn 10-50\n", +"Kc=[4 8 20]'; //different values of Kc that we have to try\n", +"K_OL=Kc*Kip*Kv*Kp*Km;//Open loop gain Eqn 10-40\n", +"K1=K_OL./(1+K_OL);//Eqn 10-38\n", +"tau1=tau./(1+K_OL); //Eqn 10-39\n", +"//========Example 11.3 now starts here========//\n", +"//Now we simulate the close loop process for these different values of K2 and tau1\n", +"M=0.05;//Magnitude of step\n", +"K2=Kp./(1+K_OL);\n", +"s=%s;\n", +"G=K2./(tau1*s+1);\n", +"G=syslin('c',G);\n", +"t=[0:0.1:10]'; //time in minutes\n", +"hdash=M*csim('step',t,G)'; \n", +"plot2d(t,hdash,rect=[0,0,10,0.2])\n", +"xgrid()\n", +"xtitle('Ex-10.3','Time(min)','h');\n", +"a=legend('Kc=4','Kc=8','Kc=20',position=4);\n", +"a.font_size=2;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;\n", +"offset=-Kp*M./(1+K_OL);\n", +"mprintf(' Kc Offset')\n", +"mprintf('\n%f %f',[Kc offset])\n", +"mprintf('\n\nNote that the book has a mistake in the question and legend of fig 10.19\n...\n", +"the values of Kc should be 4,8,20 and not 1,2,5\n...\n", +"this mistake is there because in the second edition of the book\n...\n", +"Kc has values 1 2 5 but then level transmitter span is 0.5 instead of 2')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 10.4: Stability_of_feedback_control_system.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 10.4\n", +"disp('Example 10.4')\n", +"Km=1; //We take set point gain as 1 for illustrative purposes\n", +"Kc=[15 6 2]'; //different values of Kc for which we will simulate\n", +"Gc=Kc;\n", +"s=%s;\n", +"Gv=1/(2*s+1);\n", +"Gd=1/(5*s+1);\n", +"Gp=Gd;\n", +"Gm=1/(s+1);\n", +"G=Km*Gc*Gv*Gp./(1+Km*Gc*Gv*Gp*Gm); //Eqn 10-75 G=Y/Ysp\n", +"//Now we simulate the close loop process for these different values of Kc\n", +"G=syslin('c',G);\n", +"t=[0:0.1:20]'; //time in minutes\n", +"Y=csim('step',t,G)'; \n", +"plot2d(t,Y,rect=[0,-2,20,4])\n", +"plot2d(t,ones(length(t),1),style=5)\n", +"xtitle('Ex-10.4','Time(min)','y(t)');\n", +"a=legend('Kc=15','Kc=6','Kc=2',position=3);\n", +"a.font_size=2;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/11-PID_Controller_Design_Tuning_and_Troubleshooting.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/11-PID_Controller_Design_Tuning_and_Troubleshooting.ipynb new file mode 100644 index 0000000..563760f --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/11-PID_Controller_Design_Tuning_and_Troubleshooting.ipynb @@ -0,0 +1,484 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 11: PID Controller Design Tuning and Troubleshooting" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 11.1: Direct_synthesis_for_PID.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 11.1\n", +"disp('Example 11.1')\n", +"//(a) Desired closed loop gain=1 and tau=[1 3 10]\n", +"s=%s;\n", +"tauc=[1 3 10]';\n", +"tau1=10;tau2=5;K=2;theta=1; //Time delay\n", +"Y_Ysp=(1)./(tauc*s+1); //Y/Ysp=delay/(tau*s+1) Eqn 11-6\n", +"//delay=(1-theta/2*s+theta^2/10*s^2-theta^3/120*s^3)/(1+theta/2*s+theta^2/10*s^2+theta^3/120*s^3);//Third order pade approx\n", +"delay=(1-theta/2*s)/(1+theta*s/2);//first order Pade approx\n", +"G=(K)./((tau1*s+1)*(tau2*s+1))*delay;\n", +"G_tilda=G//Model transfer function\n", +"//Eqn-11-14\n", +"Kc=1/K*(tau1+tau2)./(tauc+theta);tauI=tau1+tau2;tauD=tau1*tau2/(tau1+tau2);\n", +"Gc=Kc*(1+1/tauI/s+tauD*s); //PID without derivative filtering\n", +"G_CL=syslin('c',Gc/delay*G./(1+Gc*G));//closed loop transfer function\n", +"t=0:160;\n", +"y=csim('step',t,G_CL);\n", +"//plot(t,y)\n", +"t_d=81:160;\n", +"G_CL_dist=syslin('c',G/delay./(1+Gc*G));//closed loop wrt disturbance\n", +"u_d=[0 ones(1,length(t_d)-1)]\n", +"y_d=csim('step',t_d,G_CL_dist);\n", +"y(:,81:160)=y(:,81:160)+y_d\n", +"plot(t,y)\n", +"xgrid()\n", +"xtitle('Ex-11.1 Correct Model','Time(min)','y(t)');\n", +"a=legend('$\tau_c=1$','$\tau_c=3$','$\tau_c=10$',position=4);\n", +"a.font_size=2;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;\n", +"mprintf('\n tauc=1 tauc=3 tauc=10')\n", +"mprintf('\n Kc(K_tilda=2) %10f %f %f',Kc');\n", +"//Simulation for model with incorrect gain\n", +"scf()\n", +"K_tilda=0.9\n", +"//Eqn-11-14\n", +"Kc=1/K_tilda*(tau1+tau2)./(tauc+theta);tauI=tau1+tau2;tauD=tau1*tau2/(tau1+tau2);\n", +"Gc=Kc*(1+1/tauI/s+tauD*s) \n", +"mprintf('\n Kc(K_tilda=0.9) %10f %f %f',Kc');\n", +"mprintf('\n tauI %20f %f %f',tauI*ones(1,3));\n", +"mprintf('\n tauD %20f %f %f',tauD*ones(1,3));\n", +"G_CL=syslin('c',Gc*G./(1+Gc*G));//closed loop transfer function\n", +"t=0:160;\n", +"y=csim('step',t,G_CL);\n", +"t_d=81:160;\n", +"G_CL_dist=syslin('c',G./(1+Gc*G));//closed loop wrt disturbance\n", +"y_d=csim('step',t_d,G_CL_dist);\n", +"y(:,81:160)=y(:,81:160)+y_d\n", +"plot(t,y)\n", +"xgrid()\n", +"xtitle('Ex-11.1 Model with incorrect gain','Time(min)','y(t)');\n", +"a=legend('$\tau_c=1$','$\tau_c=3$','$\tau_c=10$',position=4);\n", +"a.font_size=2;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;\n", +"mprintf('\n \nThere is a slight mis-match between graphs from scilab code\n...\n", +"and those given in the book because of Pade approx which is very bad\n...\n", +"for delay being 1. It works only for small delays. Scilab does \n...\n", +"not handle continuous delays and hence this problem cannot \n...\n", +"be circumvented' )" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 11.3: PI_and_PID_control_of_liquid_storage_tank.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 11.3\n", +"disp('Example 11.3')\n", +"//(a)\n", +"K=0.2;theta=7.4;tauc=[8 15]';\n", +"Kc1=1/K*(2*tauc+theta)./(tauc+theta).^2; //Row M\n", +"Kc2=1/K*(2*tauc+theta)./(tauc+theta/2).^2; //Row N\n", +"tauI=2*tauc+theta;\n", +"tauD=(tauc*theta+theta^2/4)./(2*tauc+theta);\n", +"mprintf(' Kc tauI tauD')\n", +"mprintf('\nPI(tauC=8) %f %f %f',Kc1(1),tauI(1),0)\n", +"mprintf('\nPI(tauC=15) %f %f %f',Kc1(2),tauI(2),0)\n", +"mprintf('\nPID(tauC=8) %f %f %f',Kc2(1),tauI(1),tauD(1))\n", +"mprintf('\nPID(tauC=15) %f %f %f',Kc2(2),tauI(2),tauD(2))\n", +"s=%s;\n", +"//delay=(1-theta/2*s+theta^2/10*s^2-theta^3/120*s^3)/(1+theta/2*s+theta^2/10*s^2+theta^3/120*s^3);//Third order pade approx\n", +"delay=(1-theta/2*s+theta^2/10*s^2)/(1+theta/2*s+theta^2/10*s^2);//second order pade approx\n", +"//delay=(1-theta/2*s)/(1+theta/2*s);//first order pade approx\n", +"G=K*delay/s;\n", +"Gc1=Kc1.*(1+(1)./tauI/s)\n", +"Gc2=Kc2.*(1+(1)./tauI/s+tauD*s./(0.1*tauD*s+1));//PID with derivative filtering\n", +"G_CL1=syslin('c',Gc1*G./(1+Gc1*G));\n", +"G_CL2=syslin('c',Gc2*G./(1+Gc2*G));\n", +"t=0:300;\n", +"y1=csim('step',t,G_CL1);\n", +"y2=csim('step',t,G_CL2);\n", +"y1(:,1:theta)=0;//accounting for time delay--this is required otherwise\n", +"//an unrealistic inverse response is seen due to the pade approx\n", +"y2(:,1:theta)=0;\n", +"t_d=151:300;\n", +"G_CL_dist1=syslin('c',G./(1+Gc1*G));//closed loop wrt disturbance\n", +"G_CL_dist2=syslin('c',G./(1+Gc2*G));//closed loop wrt disturbance\n", +"y_d1=csim('step',t_d,G_CL_dist1);\n", +"y_d1(:,1:theta)=0;//accounting for time delay\n", +"y_d2=csim('step',t_d,G_CL_dist2);\n", +"y_d2(:,1:theta)=0;//accounting for time delay\n", +"y1(:,t_d)=y1(:,t_d)+y_d1;\n", +"y2(:,t_d)=y2(:,t_d)+y_d2;\n", +"//plot(t,y1)\n", +"//xgrid()\n", +"//xtitle('Ex-11.3 PI control','Time(min)','y(t)');\n", +"//a=legend('$\tau_c=8$','$\tau_c=15$',position=1);\n", +"//a.font_size=2;\n", +"//a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"//c=a.y_label;c.font_size=5;\n", +"//scf()\n", +"//\n", +"//plot(t,y2)\n", +"//xgrid()\n", +"//xtitle('Ex-11.3 PID control','Time(min)','y(t)');\n", +"//a=legend('$\tau_c=8$','$\tau_c=15$',position=1);\n", +"//a.font_size=2;\n", +"//a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"//c=a.y_label;c.font_size=5;\n", +"mprintf('\n\nThere is uncertainty as to whether PID with derivative filtering\n...\n", +"to be used or not. Since one gets results by using PID with filtering\n...\n", +"it has been used here. Note that pade approx for delay=7.4\n...\n", +"is totally wrong because it is too gross an approx but we have no\n...\n", +"other way of making delay approx so we have to live with it.\n\n...')\n", +"//Part (b) Routh Array testing\n", +"//For frequency response refer to ch-13 for Bode Plots\n", +"G=(1-theta*s)/s;\n", +"poly_PI=Gc1*G;//denom(G_CL1);//G*Gc for PI controller\n", +"poly_PID=Gc2*G;//G*Gc for PID controller\n", +"Routh1=routh_t(poly_PI(1,1)/1,poly(0,'K')); // produces routh table for polynomial 1+Kc*poly\n", +"disp(Routh1,'Routh1=')\n", +"Kmax1=roots(numer(Routh1(1,1)));\n", +"Routh2=routh_t(poly_PI(2,1)/1,poly(0,'K')); // produces routh table for polynomial 1+Kc*poly\n", +"disp(Routh2,'Routh2=')\n", +"Kmax2=roots(numer(Routh2(1,1)));\n", +"Routh3=routh_t(poly_PID(1,1)/1,poly(0,'K')); // produces routh table for polynomial 1+Kc*poly\n", +"disp(Routh3,'Routh3=')\n", +"//Kmax3=roots(numer(Routh3(1,1)));\n", +"Routh4=routh_t(poly_PID(2,1)/1,poly(0,'K')); // produces routh table for polynomial 1+Kc*poly\n", +"disp(Routh4,'Routh4=')\n", +"//Kmax4=roots(numer(Routh4(1,1)));\n", +"mprintf('\n Kmax should be less than %f and %f \n for tauc=8 and 15 respectively for PI system to be stable',Kmax1,Kmax2)\n", +"mprintf('\n\nAnswers to Kmax for PID controller using \n...\n", +"Routh Array in the book are wrong. This can be easily \n...\n", +"checked from Routh3 and Routh4 which are displayed\n')\n", +"mprintf('\n\nFor frequency response refer to ch-13 for Bode Plots\n')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 11.4: IMC_for_lag_dominant_model.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 11.4\n", +"disp('Example 11.4')\n", +"s=%s;\n", +"theta=1;tau=100;K=100;\n", +"delay=(1-theta/2*s+theta^2/10*s^2-theta^3/120*s^3)/(1+theta/2*s+theta^2/10*s^2+theta^3/120*s^3);//Third order pade approx\n", +"G=K*delay/(tau*s+1);\n", +"//(a)\n", +"tauca=1;\n", +"Kc1=1/K*tau/(tauca+theta);taui1=tau;\n", +"//(b)\n", +"taucb=2;Kstar=K/tau;\n", +"Kc2=1/Kstar*(2*taucb+theta)./(taucb+theta).^2; //Row M\n", +"taui2=2*taucb+theta;\n", +"//(c)\n", +"taucc=1;\n", +"Kc3=Kc1;taui3=min(taui1,4*(taucc+theta))\n", +"//(d)\n", +"//Kc4=0.551;taui4=4.91;\n", +"//Chen and Seborg settings given in Second Edition of book\n", +"mprintf(' Kc tauI')\n", +"mprintf('\nIMC %20f %f',Kc1,taui1)\n", +"mprintf('\nIntegrator approx %-5f %f',Kc2,taui2)\n", +"mprintf('\nSkogestad %15f %f',Kc3,taui3)\n", +"//mprintf('\nDS-d %20f %f',Kc4,taui4) \n", +"Gc=[Kc1 Kc2 Kc3]'.*(1+(1)./([taui1 taui2 taui3]'*s))\n", +"G_CL=syslin('c',Gc*G./(1+Gc*G));\n", +"t=0:0.1:20;\n", +"y=csim('step',t,G_CL);\n", +"y(:,1:theta/0.1)=0;//accounting for time delay--this is required otherwise\n", +"//an unrealistic inverse response is seen due to the pade/taylor approx\n", +"plot(t,y);\n", +"xgrid()\n", +"xtitle('Ex-11.4 Tracking problem','Time(min)','y(t)');\n", +"a=legend('IMC','Integrator approx','Skogestad',position=4);\n", +"a.font_size=2;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;\n", +"scf()\n", +"t=0:0.1:60;\n", +"G_CL_dist=syslin('c',G./(1+Gc*G));//closed loop wrt disturbance\n", +"yd=csim('step',t,G_CL_dist);\n", +"yd(:,1:theta/0.1)=0;//accounting for time delay--this is required otherwise\n", +"//an unrealistic inverse response is seen due to the pade/taylor approx\n", +"plot(t,yd);\n", +"xgrid()\n", +"xtitle('Ex-11.4 Disturbance rejection','Time(min)','y(t)');\n", +"a=legend('IMC','Integrator approx','Skogestad',position=4);\n", +"a.font_size=2;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 11.5: PI_controller_IMC_ITAE.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 11.5\n", +"disp('Example 11.5')\n", +"K=1.54;theta=1.07;tau=5.93;\n", +"//(a)\n", +"tauca=tau/3;\n", +"Kc1=1/K*tau/(tauca+theta);taui1=tau; //Table 11.1\n", +"//(b)\n", +"taucb=theta;\n", +"Kc2=1/K*tau/(taucb+theta);taui2=tau; //Table 11.1\n", +"//(c)\n", +"//Table 11.3\n", +"Y=0.859*(theta/tau)^(-0.977);Kc3=Y/K;\n", +"taui3=tau*inv(0.674*(theta/tau)^-0.680); \n", +"//(d)\n", +"//Table 11.3\n", +"Kc4=1/K*0.586*(theta/tau)^-0.916;taui4=tau*inv(1.03-0.165*(theta/tau)); \n", +"mprintf(' Kc tauI')\n", +"mprintf('\nIMC(tauC=tau/3) %f %f',Kc1,taui1)\n", +"mprintf('\nIMC(tauC=theta) %-5f %f',Kc2,taui2)\n", +"mprintf('\nITAE(disturbance) %f %f',Kc3,taui3)\n", +"mprintf('\nITAE(set point) %10f %f',Kc4,taui4)" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 11.6: Controller_with_two_degrees_of_freedom.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 11.6\n", +"disp('Example 11.6')\n", +"//Drawing on example 11.4\n", +"s=%s;\n", +"theta=1;tau=100;K=100;\n", +"delay=(1-theta/2*s+theta^2/10*s^2-theta^3/120*s^3)/(1+theta/2*s+theta^2/10*s^2+theta^3/120*s^3);//Third order pade approx\n", +"G=K*delay/(tau*s+1);\n", +"Kc=0.556;taui=5;\n", +"Gc=Kc.*(1+(1)./([taui]*s))\n", +"G_CL=syslin('c',Gc*G./(1+Gc*G));\n", +"t=0:0.1:20;\n", +"y1=csim('step',t,G_CL);\n", +"y1(:,1:theta/0.1)=0;//accounting for time delay--this is required otherwise\n", +"//an unrealistic inverse response is seen due to the pade/taylor approx\n", +"beta=0.5;\n", +"G_CL2=syslin('c',(Gc+beta-1)*G./(1+Gc*G));//This can be obtained on taking\n", +"//laplace transform of eqn 11-39 and making a block diagram \n", +"//In Eqn 11-39 p refers to input to the process\n", +"t=0:0.1:20;\n", +"y2=csim('step',t,G_CL2);\n", +"y2(:,1:theta/0.1)=0;//accounting for time delay--this is required otherwise\n", +"//an unrealistic inverse response is seen due to the pade/taylor approx\n", +"plot(t,[y1; y2]);\n", +"xgrid()\n", +"xtitle('Ex-11.3 Tracking problem','Time(min)','y(t)');\n", +"a=legend('$\beta=1$','$\beta=0.5$',position=4);\n", +"a.font_size=2;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;\n", +"//Note that there is a slight mis-match between the plots obtained from scilab code\n", +"//and that of the book because of third order pade approximation\n", +"//The plots in the book have been produced using advanced proprietary software\n", +"//which supports using exact delays while scilab does not have that functionality" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 11.7: Continuous_cycling_method.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 11.7\n", +"disp('Example 11.7')\n", +"K=2;theta=1;tau1=10;tau2=5; //Model parameters\n", +"s=%s;\n", +"delay=(1-theta/2*s+theta^2/10*s^2-theta^3/120*s^3)/(1+theta/2*s+theta^2/10*s^2+theta^3/120*s^3);//Third order pade approx\n", +"G=K*delay/((tau1*s+1)*(tau2*s+1));\n", +"Ku=[8.01]';//Trials for various values of Ku can be done by changing Ku\n", +"G_CL_trial=syslin('c',G*Ku./(1+G*Ku))\n", +"t=0:0.1:100;\n", +"y=csim('step',t,G_CL_trial);\n", +"plot(t,y)\n", +"xtitle('Ex-11.7 Finding ultimate gain Ku','Time(min)','y(t)');\n", +"a=legend('Closed loop test',position=4);\n", +"a.font_size=2;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;\n", +"//There isnot a sustained oscillation for Ku=7.88, in our simulation because\n", +"//we are using a third order Pade Approx for delay. But still we go ahead with it\n", +"//so that it matches with the example values. Our simulations give Ku=8\n", +"Ku=7.88;Pu=11.66;\n", +"//(a) Table 11.4 ZN\n", +"Kc1=0.6*Ku;taui1=Pu/2;tauD1=Pu/8;\n", +"//(b) Table 11.4 TL\n", +"Kc2=0.45*Ku;taui2=Pu*2.2;tauD2=Pu/6.3;\n", +"//(c) DS method\n", +"tauc=3;\n", +"Kc3=1/K*(tau1+tau2)/(tauc+theta);taui3=tau1+tau2;tauD3=tau1*tau2/(tau1+tau2);\n", +"mprintf(' Kc tauI tauD')\n", +"mprintf('\nZN %f %f %f',Kc1,taui1,tauD1)\n", +"mprintf('\nTL %f %f %f',Kc2,taui2,tauD2)\n", +"mprintf('\nDS %f %f %f',Kc3,taui3,tauD3)\n", +"// Now we compare the performance of each controller setting\n", +"Gc1=Kc1.*(1+(1)./taui1/s+tauD1*s)\n", +"Gc2=Kc2.*(1+(1)./taui2/s+tauD2*s)\n", +"Gc3=Kc3.*(1+(1)./taui3/s+tauD3*s)\n", +"Gc=[Gc1 Gc2 Gc3]';\n", +"G_CL=syslin('c',Gc*G./(1+Gc*G));\n", +"t=0:160;\n", +"y=csim('step',t,G_CL);\n", +"y(:,1:theta)=0;//accounting for time delay--this is required otherwise\n", +"//an unrealistic inverse response is seen due to the pade/taylor approx\n", +"t_d=81:160;\n", +"G_CL_dist=syslin('c',G./(1+Gc*G));//closed loop wrt disturbance\n", +"yd=csim('step',t_d,G_CL_dist);\n", +"yd(:,1:theta)=0;//accounting for time delay\n", +"y(:,t_d)=y(:,t_d)+yd;\n", +"scf();\n", +"plot(t,y)\n", +"xgrid()\n", +"xtitle('Ex-11.7 Comparison of 3 controllers','Time(min)','y(t)');\n", +"a=legend('Ziegler-Nichols','Tyerus-Luyben','Direct Synthesis',position=4);\n", +"a.font_size=2;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;\n", +"mprintf('\n\nThere is slight mismatch between scilab simulation\n...\n", +"and book simulation due to Pade approx\n')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 11.8: Reaction_curve_method.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 11.8\n", +"disp('Example 11.8')\n", +"S=(55-35)/(7-1.07);//%/min\n", +"delta_p=43-30;//%\n", +"R=S/delta_p;//min^-1\n", +"delta_x=55-35;//%\n", +"K=delta_x/delta_p;\n", +"theta=1.07;//min\n", +"tau=7-theta;//min\n", +"mprintf('\nThe resulting process model is with delay of 1.07 min\n')\n", +"s=%s;\n", +"G=K/(tau*s+1);\n", +"disp(G,'G=')" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/12-Control_Strategies_at_the_Process_Unit_Level.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/12-Control_Strategies_at_the_Process_Unit_Level.ipynb new file mode 100644 index 0000000..1996cfc --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/12-Control_Strategies_at_the_Process_Unit_Level.ipynb @@ -0,0 +1,59 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 12: Control Strategies at the Process Unit Level" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 12.1: Degrees_of_freedom.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 12.1\n", +"disp('Example 12.1')\n", +"NE=3;\n", +"NV=6;\n", +"NF=NV-NE;\n", +"ND=2;\n", +"NFC=NF-ND;\n", +"mprintf(' NF=%i\n NFC=%i',NF,NFC)" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/13-Frequency_response_analysis_and_control_system_design.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/13-Frequency_response_analysis_and_control_system_design.ipynb new file mode 100644 index 0000000..e8ecc19 --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/13-Frequency_response_analysis_and_control_system_design.ipynb @@ -0,0 +1,433 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 13: Frequency response analysis and control system design" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 13.3: Bode_Plot.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 13.3\n", +"disp('Example 13.3')\n", +"function bodegen(num,den,w,lf,delay)\n", +"//Bode plot\n", +"//Numerator and denominator are passed as input arguments\n", +"//Both are polynomials in powers of s(say)\n", +"//This function has been modified from the original one \n", +"//written by Prof Kannan Moudgalya, Dept of ChemE, IIT-B\n", +"G = num/den;\n", +"G1 = horner(G,%i*w);\n", +"G1p = phasemag(G1)-delay*w*180/%pi;\n", +"if LF == 'normal' then \n", +" xset('window',0); clf();\n", +" subplot(2,1,1)\n", +" plot2d(w,abs(G1),logflag='nn',style = 2); \n", +" xtitle('Normal scale','','Magnitude'); xgrid();\n", +" subplot(2,1,2)\n", +" plot2d1(w,G1p,logflag='nn',style = 2);\n", +" xgrid();\n", +" xtitle('w(rad/sec)','','Phase(deg)');\n", +"elseif LF == 'semilog' then \n", +" xset('window',1); clf();\n", +" subplot(2,1,1)\n", +" plot2d(w,20*log10(abs(G1)),logflag='ln',style = 2);\n", +" xgrid(); \n", +" xtitle('Semilog','','Magnitude (dB)');\n", +" subplot(2,1,2)\n", +" plot2d1(w,G1p,logflag='ln',style = 2);\n", +" xgrid();\n", +" xtitle('w(rad/sec)','','Phase(deg)');\n", +"elseif LF == 'loglog' then \n", +" xset('window',2); clf();\n", +" subplot(2,1,1)\n", +" plot2d(w,abs(G1),logflag='ll',style = 2); \n", +" xgrid();\n", +" xtitle('Loglog','','Magnitude');\n", +" subplot(2,1,2)\n", +" plot2d1(w,G1p,logflag='ln',style = 2,rect=[0.01,-270,100,0]);//note the usage of rect for this particular example\n", +" xgrid();\n", +" xtitle('w(rad/sec)','','Phase(deg)');\n", +"end\n", +"endfunction;\n", +"s = %s;\n", +"num = 5*(0.5*s+1);\n", +"den = (20*s+1)*(4*s+1);\n", +"theta=1;\n", +"w = 0.001:0.002:10*%pi;\n", +"LF = 'loglog' // Warning: Change this as necessary\n", +"bodegen(num,den,w,LF,theta);\n", +"//Checking using iodelay toolbox in scilab\n", +"//G=syslin('c',num/den);\n", +"//G=iodelay(G,1)\n", +"//bode(G,0.01,0.1)" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 13.4: Bode.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 13.4\n", +"disp('Example 13.4')\n", +"s = %s;\n", +"num = 2;\n", +"den = (0.5*s+1)^3;\n", +"delay=0;\n", +"w = 0.001:0.002:100*%pi;\n", +"LF = 'loglog' // Warning: Change this as necessary\n", +"//Kc=1\n", +"G1 = num/den;\n", +"G1m = horner(G1,%i*w); //G1m denotes magnitude\n", +"G1p = phasemag(G1m)-delay*w*180/%pi; //G1p denotes phase\n", +"//Kc=4\n", +"G2 = 4*num/den;\n", +"G2m = horner(G2,%i*w);\n", +"G2p = phasemag(G2m)-delay*w*180/%pi;\n", +"//Kc=20\n", +"G3 = 20*num/den;\n", +"G3m = horner(G3,%i*w);\n", +"G3p = phasemag(G3m)-delay*w*180/%pi;\n", +" xset('window',0); \n", +" subplot(2,1,1)\n", +" plot2d(w,abs(G3m),logflag='ll',style = 4);\n", +" plot2d(w,abs(G2m),logflag='ll',style = 3);\n", +" plot2d(w,abs(G1m),logflag='ll',style = 2); \n", +" \n", +" xgrid();\n", +" xtitle('Loglog','','Magnitude'); \n", +" legend('Kc=20','Kc=4','Kc=1')\n", +" subplot(2,1,2)\n", +" plot2d1(w,G1p,logflag='ln',style = 2);\n", +" xgrid();\n", +" xtitle('w(rad/sec)','','Phase(deg)');" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 13.5: PI_control_of_overdamped_second_order_process.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 14.5\n", +"disp('Example 14.5')\n", +"s = %s;\n", +"num = 1;\n", +"den = (9*s+1)*(11*s+1);\n", +"delay=0.3;\n", +"w = 0.001:0.002:5*%pi;\n", +"LF = 'loglog' // Warning: Change this as necessary\n", +"Gc=20*(1+1/2.5/s+s);\n", +"G1 = num/den*Gc;\n", +"G1m = horner(G1,%i*w); //G1m denotes magnitude\n", +"G1p = phasemag(G1m)-delay*w*180/%pi; //G1p denotes phase\n", +" xset('window',0); \n", +" subplot(2,1,1)\n", +" plot2d(w,abs(G1m),logflag='ll',style = 3);\n", +" xgrid();\n", +" xtitle('Loglog','','Magnitude'); \n", +" subplot(2,1,2)\n", +" plot2d1(w,G1p,logflag='ln',style = 1);\n", +" xgrid();\n", +" xtitle('w(rad/sec)','','Phase(deg)');" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 13.6: Bode_plot.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 13.6\n", +"disp('Example 13.6')\n", +"s = %s;\n", +"num = 1;\n", +"den = (9*s+1)*(11*s+1);\n", +"delay=0.3;\n", +"w = 0.001:0.002:5*%pi;\n", +"LF = 'loglog' // Warning: Change this as necessary\n", +"Gc=20*(1+1/2.5/s+s);\n", +"G1 = num/den*Gc;\n", +"G1m = horner(G1,%i*w); //G1m denotes magnitude\n", +"G1p = phasemag(G1m)-delay*w*180/%pi; //G1p denotes phase\n", +" xset('window',0); \n", +" subplot(2,1,1)\n", +" plot2d(w,abs(G1m),logflag='ll',style = 3);\n", +" xgrid();\n", +" xtitle('Loglog','','Magnitude'); \n", +" subplot(2,1,2)\n", +" plot2d1(w,G1p,logflag='ln',style = 1);\n", +" xgrid();\n", +" xtitle('w(rad/sec)','','Phase(deg)');" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 13.7: Bode_Plot.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 13.7\n", +"disp('Example 13.7')\n", +"s = %s;\n", +"num = 4;\n", +"den = (5*s+1);\n", +"delay=1;\n", +"w = 0.001:0.002:10*%pi;\n", +"LF = 'loglog' // Warning: Change this as necessary\n", +"Gv=2;Gm=0.25;Gc=1;\n", +"G1 = num/den*Gc*Gm*Gv;\n", +"G1m = horner(G1,%i*w); //G1m denotes magnitude\n", +"G1p = phasemag(G1m)-delay*w*180/%pi; //G1p denotes phase\n", +" xset('window',0); \n", +" subplot(2,1,1)\n", +" plot2d(w,abs(G1m),logflag='ll',style = 3,rect=[0.01,0.1,100,10]);\n", +" xgrid();\n", +" xtitle('Loglog','','Magnitude'); \n", +" subplot(2,1,2)\n", +" plot2d1(w,G1p,logflag='ln',style = 1,rect=[0.01,-400,100,100]);\n", +" xgrid();\n", +" xtitle('w(rad/sec)','','Phase(deg)');\n", +" \n", +"//Example ends\n", +" \n", +"//bin/\n", +"//boot/\n", +"//dev/\n", +"//etc/\n", +"//home/\n", +"//lib/\n", +"//lib32/\n", +"//lib64/\n", +"//libx32/\n", +"//lost+found/\n", +"//media/\n", +"//mnt/\n", +"//opt/\n", +"//proc/\n", +"//root/\n", +"//run/\n", +"//sbin/\n", +"//srv/\n", +"//sys/\n", +"//tmp/\n", +"//usr/\n", +"//var/\n", +"//In the SECOND EDITION of the book, this example also asks for drawing Nyquist plot\n", +"//In case you want to learn how to do it, Uncomment the code below\n", +" \n", +"////Please install IODELAY toolbox from Modeling and Control tools in ATOMS\n", +"////http://atoms.scilab.org/toolboxes/iodelay/0.4.5\n", +"////There is no inbuilt toolbox in scilab for introducing time delays other than\n", +"////above mentioned. The output of iodelay toolbox however does not work\n", +"////with csim and syslin commands\n", +"////The output of iodelay however can be used for frequency related analyses\n", +"////like bode and nyquist\n", +"//xset('window',1);\n", +"//G=syslin('c',G1);\n", +"//G=iodelay(G,delay);\n", +"//nyquist(G,%f); //%f => asymmetric, see help nyquist\n", +"////bin/\n", +"////boot/\n", +"////dev/\n", +"////etc/\n", +"////home/\n", +"////lib/\n", +"////lib32/\n", +"////lib64/\n", +"////libx32/\n", +"////lost+found/\n", +"////media/\n", +"////mnt/\n", +"////opt/\n", +"////proc/\n", +"////root/\n", +"////run/\n", +"////sbin/\n", +"////srv/\n", +"////sys/\n", +"////tmp/\n", +"////usr/\n", +"////var/" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 13.8: Maximum_permissible_time_delay_for_stability.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 13.8\n", +"disp('Example 13.8')\n", +"s = %s;\n", +"num = 4;\n", +"den = (5*s+1);\n", +"delay=1;\n", +"w = 0.001:0.002:10*%pi;\n", +"LF = 'loglog' // Warning: Change this as necessary\n", +"Gv=2;Gm=0.25;\n", +"Ku=4.25;Pu=2*%pi/1.69;\n", +"//Ziegler Nichols\n", +"Kc1=0.6*Ku;taui1=Pu/2;tauD1=Pu/8;\n", +"//Tyreus Luyben\n", +"Kc2=0.45*Ku;taui2=Pu*2.2;tauD2=Pu/6.3;\n", +"mprintf(' Kc tauI tauD')\n", +"mprintf('\nZN %f %f %f',Kc1,taui1,tauD1)\n", +"mprintf('\nTL %f %f %f',Kc2,taui2,tauD2)\n", +"Gc_ZN=Kc1*(1+1/taui1/s+s*tauD1/(0.1*s*tauD1+1));\n", +"Gc_TL=Kc2*(1+1/taui2/s+s*tauD2/(0.1*s*tauD2+1)); //Filtered Controllers with filter constant as 0.1\n", +"G1 = num/den*Gc_ZN*Gm*Gv;\n", +"G1m = horner(G1,%i*w); //G1m denotes magnitude\n", +"Abs_G1m=abs(G1m)\n", +"G1p = phasemag(G1m)-delay*w*180/%pi; //G1p denotes phase\n", +"G2 = num/den*Gc_TL*Gm*Gv;\n", +"G2m = horner(G2,%i*w); //G1m denotes magnitude\n", +"Abs_G2m=abs(G2m);\n", +"G2p = phasemag(G2m)-delay*w*180/%pi; //G1p denotes phase\n", +" xset('window',0); \n", +" subplot(2,1,1)\n", +" plot2d(w,Abs_G1m,logflag='ll',style = 1,rect=[0.01,0.1,10,1000]);\n", +" plot2d(w,Abs_G2m,logflag='ll',style = 2,rect=[0.01,0.1,10,1000]);\n", +" legend('Ziegler-Nichols','Tyreus Luyben')\n", +" xgrid();\n", +" xtitle('Loglog','','Magnitude'); \n", +" subplot(2,1,2)\n", +" plot2d1(w,G1p,logflag='ln',style = 1,rect=[0.01,-300,10,-50]);\n", +" plot2d1(w,G2p,logflag='ln',style = 2,rect=[0.01,-300,10,-50]);\n", +"// legend('Ziegler-Nichols','Tyreus Luyben')\n", +" xgrid();\n", +" xtitle('w(rad/sec)','','Phase(deg)');\n", +" \n", +"//G_ZN=syslin('c',G1);\n", +"//G_ZN=iodelay(G_ZN,delay);\n", +"//G_TS=syslin('c',G2);\n", +"//G_TS=iodelay(G_TS,delay);\n", +"//scf();nyquist(G_TS,%f)\n", +"//[gm_ZN,fr_ZN]=g_margin(G_ZN);[gm_TS,fr_TS]=g_margin(G_TS);\n", +"//[pm_ZN,fr_ZN_p]=p_margin(G_ZN);[pm_ZN,fr_ZN_p]=p_margin(G_TS);\n", +"//g_maring and p_margin do not support iodelay toolbox, hence we \n", +"//cannot use these so we try a workaround\n", +"//We can find w for which magnitude(AR) is 1 and \n", +"//and calculate phase corresponding to it which gives us phase margin\n", +"//Also we can find crossover frequency and thus find Gain Margin\n", +" indices1=find(abs(Abs_G1m-1)<0.01) //We find those values of indices of Abs_G1m for which it is almost 1\n", +" indices2=find(abs(Abs_G2m-1)<0.01)\n", +" //size(indices)\n", +" PM1=mean(G1p(indices1))+180\n", +" PM2=mean(G2p(indices2))+180\n", +" \n", +" indices1_p=find(abs(G1p+180)<0.05) //We find those values of indices of G1p for which it is almost -180\n", +" indices2_p=find(abs(G2p+180)<0.05)\n", +" //size(indices)\n", +" GM1=1/mean(Abs_G1m(indices1_p))\n", +" GM2=1/mean(Abs_G2m(indices2_p))\n", +" \n", +" wc1=mean(w(indices1_p)); \n", +" wc2=mean(w(indices2_p));\n", +"mprintf('\n\n\n GM PM wc')\n", +"mprintf('\nZN %f %f %f',GM1,PM1,wc1)\n", +"mprintf('\nTL %f %f %f\n',GM2,PM2,wc2)\n", +"theta=PM2*%pi/0.79/180;\n", +"disp(theta,'deltatheta(min)=')" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/14-Feedforward_and_Ratio_Control.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/14-Feedforward_and_Ratio_Control.ipynb new file mode 100644 index 0000000..b1e0887 --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/14-Feedforward_and_Ratio_Control.ipynb @@ -0,0 +1,182 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 14: Feedforward and Ratio Control" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 14.1: Ratio_control.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 14.1\n", +"disp('Example 14.1')\n", +"Sd=30;\n", +"Su=15;\n", +"Rd=1/3;\n", +"K_R=Rd*Sd/Su; //Eqn 14-3\n", +"mprintf(' K_R=%f',K_R)" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 14.5: Feedforward_control_in_blending_process.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 14.5\n", +"disp('Example 14.5')\n", +"mprintf('\n\nThere are many errors in this example\n...\n", +"1.In Eqn 14-17 the value of w2_o is not zero. It is 50kg/min.\n...\n", +"This is so because otherwise current signal from p(t) ie ...\n", +"\n eqn 14-30 is more than 20mA which is not possible\n\n....\n", +"2.The step change in x1 is from 0.2 to 0.3 and not 0.2 to 0.4\n...\n", +"If there is a step change to x1=0.4, then with x2=0.6\n...\n", +"one cannot achieve output xsp=0.34 because it is less\n...\n", +"both x1 and x2.\n\n...\n", +"3.The gain of Gd is 0.65 which is correct because V\n...\n", +"has to be calculated using height=1.5meter ie\n...\n", +"how much the CSTR is filled and not h=3m, ie\n...\n", +"the capacity of CSTR. This is important because\n...\n", +"the person who has made solutions for the book has taken h=3m\n...\n", +"for generating graphs and hence the gain is twice. \n...\n", +"the graphs generated from this code are correct\n\n' )\n", +"//part(a) //========Static feedforward controller==========//\n", +"K_IP=(15-3)/(20-4);\n", +"Kv=300/12;tauV=0.0833;\n", +"Kt=(20-4)/0.5;\n", +"w2_o=50;x1_o=0;//Zero of the instrument\n", +"w1bar=650;w2bar=350;//kg/min\n", +"C1=4-w2_o/Kv/K_IP; //Eqn 14-16 to 14-19\n", +"C2=w1bar/(Kv*K_IP*Kt);\n", +"C3=4+Kt*x1_o;\n", +"x1bar=0.2;x2bar=0.6;xbar=0.34;\n", +"mprintf('\nThe values of C1, C2, C3 in Eqns 14-16 to 14-19 are\n %f, %f, %f',C1,C2,C3)\n", +"//part(b) //=======Dynamic feedforward controller=======//\n", +"s=%s;\n", +"theta=1;\n", +"V=%pi*1^2*1.5; //pi*r^2*h finding volume\n", +"rho=1000;//kg/m3\n", +"wbar=w1bar+w2bar;\n", +"tauD=V*rho/w2bar;tauP=V*rho/wbar;\n", +"Kp=(x2bar-xbar)/wbar;\n", +"Kd=w1bar/wbar;\n", +"Gv=Kv/(tauV*s+1);\n", +"Gd=Kd/(tauP*s+1);\n", +"Gt=Kt;delay=1;\n", +"Gp=Kp/(tauP*s+1);\n", +"Gf=-Gd/Gv/Gt/Gp/K_IP; //Equation 14-33 without exp(+s)\n", +"//Gt=32*(1-theta/2*s+theta^2/12*s^2)/(1+theta/2*s+theta^2/12*s^2);//second order Pade approx.\n", +"Gt=32*(1-theta/2*s)/(1+theta/2*s);//first order Pade approx.\n", +"alpha=0.1;\n", +"Gf=horner(Gf,0)*(1.0833*s+1)/(alpha*1.0833*s+1);//Eqn 14-34\n", +"disp(Gf,'Gf=')\n", +"//========Static feedforward controller simulation==========//\n", +"Ts=0.01;//sampling time in minutes\n", +"t=Ts:Ts:50;\n", +"xsp=0.34;//set point for conc. output of blender\n", +"x1step=0.2+[zeros(1,length(t)/5) 0.1*ones(1,length(t)*4/5)];//disturbance\n", +"//There is a one second lag in the measurement of the disturbance by Gt\n", +"delay=1;\n", +"d=[0.2*ones(1,delay/Ts) x1step(1,1:$-delay/Ts)];//measurement of disturbance with lag\n", +"x1m=4+Kt*d; //Eqn 14-12 where d=x1(t)-x1_o\n", +"//plot(t,[x1step' x1m'])\n", +"pt=C1+C2*(Kt*xsp-x1m+C3)/(x2bar-xsp);\n", +"//Now the values calculated by the above controller needs to be passed to the process \n", +"G_static=syslin('c',[Gd K_IP*Gv*Gp]);\n", +"//we take disturbance and control action in deviation variables\n", +"yt=0.34+csim([x1step-x1step(1,1);pt-pt(1,1)],t,G_static);\n", +"subplot(2,1,1)\n", +"plot(t,yt);\n", +"xtitle('Fig 14.13(a)','Time(min)','x(mass fraction)')\n", +"xgrid();\n", +"//========Dynamic feedforward controller simulation======//\n", +"Ys_Ds=[Gd K_IP*32*Gf*Gv*Gp]; //Gt=32 without delay\n", +"Ys_Ds=syslin('c',Ys_Ds);\n", +"t=0.01:0.01:50;\n", +"d=[zeros(1,length(t)/5) 0.1*ones(1,length(t)*4/5)];//disturbance\n", +"d_shift=[zeros(1,1.1*length(t)/5) 0.1*ones(1,length(t)*3.9/5)];\n", +"//we delay the disturbance by one minute for the feedforward controller\n", +"//We do this because Pade approx is not good for delay of 1 minute\n", +"yt=0.34+csim([d;d_shift],t,Ys_Ds)\n", +"plot(t,yt,color='red')\n", +"legend('Static Gain','Dynamic Compensation')\n", +"//part(c) //=======PI controller for Feedback======//\n", +"Kcu=48.7;Pu=4;//min\n", +"Kc=0.45*Kcu;tauI=Pu/1.2;tauD=0;\n", +"Gc=Kc*(1+1/(tauI*s)+tauD*s/(1+tauD*s*0.1));\n", +"Gm=Gt;//sensor/transmitter\n", +"//==========Feedforward and feedback control with dynamic compensation======//\n", +"Ys_Ds=[Gd K_IP*32*Gf*Gv*Gp]/(1+K_IP*Gc*Gv*Gp*Gm);//32 is magnitude of Gt\n", +"Ys_Ds=syslin('c',Ys_Ds);\n", +"t=0.01:0.01:50;\n", +"d=[zeros(1,length(t)/5) 0.1*ones(1,length(t)*4/5)];//disturbance\n", +"yt=0.34+csim([d;d_shift],t,Ys_Ds)\n", +"//This shifting is better because Pade approx is not accurate. Note that there is \n", +"//pade approx in the denominator also(Gm) which we cant help.\n", +"subplot(2,1,2)\n", +"plot(t,yt)\n", +"xgrid();\n", +"xtitle('Fig 14.13(b)','Time(min)','x(mass fraction)')\n", +"//==========Feedback control only with dynamic compensation======//\n", +"Ys_Ds=(Gd)/(1+K_IP*Gc*Gv*Gp*Gm);\n", +"Ys_Ds=syslin('c',Ys_Ds);\n", +"d=[zeros(1,length(t)/5) 0.1*ones(1,length(t)*4/5)];//disturbance\n", +"yt=0.34+csim(d,t,Ys_Ds)\n", +"plot(t,yt,color='red')\n", +"legend('FB+FF with dynamic compensation','FB only')\n", +"mprintf('\n\nNote the slight mismatch between response \n...\n", +"times due to pade approx the gain is half of that in the\n...\n", +"book. Please see the heigh explanation above to understand.')" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/15-Enhanced_Single_Loop_Control_Strategies.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/15-Enhanced_Single_Loop_Control_Strategies.ipynb new file mode 100644 index 0000000..9ba3aa2 --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/15-Enhanced_Single_Loop_Control_Strategies.ipynb @@ -0,0 +1,114 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 15: Enhanced Single Loop Control Strategies" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 15.1: Stability_limits_for_proportional_cascade_controller.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 15.1\n", +"disp('Example 15.1')\n", +"s=%s;\n", +"Gp1=4/((4*s+1)*(2*s+1));Gp2=1;Gd2=1;Gd1=1/(3*s+1);\n", +"Gm1=0.05;Gm2=0.2;\n", +"Gv=5/(s+1);\n", +"Kc2=4;\n", +"Ys=Kc2*Gv*Gp1*Gm1/(1+Kc2*Gv*Gm2);\n", +"Routh=routh_t(Ys,poly(0,'Kc1')); // produces routh table for polynomial 1+Kc*Ys\n", +"disp(Routh)\n", +"K1=roots(numer(Routh(3,1)));\n", +"K2=roots(numer(Routh(4,1)));\n", +"mprintf('\n Kc1 lies between %f and %f \n for cascade system to be stable', K2,K1)\n", +"Ys2=Gv*Gp2*Gp1*Gm1;\n", +"Routh2=routh_t(Ys2,poly(0,'Kc1')); // produces routh table for polynomial 1+Kc*Ys\n", +"disp(Routh2)\n", +"K1_2=roots(numer(Routh2(3,1)));\n", +"K2_2=roots(numer(Routh2(4,1)));\n", +"mprintf('\n Kc1 lies between %f and %f \n for conventional system to be stable', K2_2,K1_2)\n", +"//We cannot find offset symbolically in Scilab because scilab does not support\n", +"//handling of two variables in single polynomial\n", +"//To find this limit one can use Sage\n", +"//However in this case the calculations can be done in a very easy way by hand\n", +"//and hence do not require to be computed from Sage" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 15.2: Set_point_response_for_second_order_transfer_function.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 15.2\n", +"disp('Example 15.2')\n", +"s = %s;\n", +"theta=1 //delay\n", +"delay=(1-theta/2*s+theta^2/12*s^2)/(1+theta/2*s+theta^2/12*s^2);//Second order pade approx\n", +"G=1/((5*s+1)*(3*s+1));\n", +"Gp=[G;delay*G];//Both models with and without delay\n", +"Gc=[3.02*(1+1/(6.5*s));1.23*(1+1/(7*s))];\n", +"G_CL=syslin('c',(Gp.*Gc)./(1+Gp.*Gc))\n", +"t=0:0.01:40;\n", +"yt=csim('step',t,G_CL)\n", +"plot2d(t',yt') //For plotting multiple graphs in one command make sure time is n*1 vector \n", +"//while yt is n*p vector where p are the no. of plots\n", +"xtitle('Example-15.2','Time(min)','$y(t)$');\n", +"xgrid(); \n", +"a=legend('delay=0','delay=1',position=4);\n", +"a.font_size=2;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/16-Multiloop_and_Multivariable_Control.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/16-Multiloop_and_Multivariable_Control.ipynb new file mode 100644 index 0000000..6755242 --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/16-Multiloop_and_Multivariable_Control.ipynb @@ -0,0 +1,335 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 16: Multiloop and Multivariable Control" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 16.1: Pilot_scale_distillation_column.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 16.1\n", +"disp('Example 16.1')\n", +"K=[12.8 -18.9;6.6 -19.4];\n", +"tau=[16.7 21;10.9 14.4];\n", +"s=%s;\n", +"G=K./(1+tau*s);\n", +"//ITAE settings from Table 11.3\n", +"K1=12.8;tau1=16.7;theta1=1;K2=-19.4;tau2=14.4;theta2=3;\n", +"Kc1=1/K1*0.586*(theta1/tau1)^-0.916;taui1=tau1*inv(1.03-0.165*(theta1/tau1)); \n", +"Kc2=1/K2*0.586*(theta2/tau2)^-0.916;taui2=tau2*inv(1.03-0.165*(theta2/tau2)); \n", +"mprintf(' Kc tauI')\n", +"mprintf('\nx_D-R %f %f',Kc1,taui1)\n", +"mprintf('\nx_B-R %f %f',Kc2,taui2)\n", +"Kc=[Kc1;Kc2];\n", +"tauI=[taui1;taui2];\n", +"//====Making step response models of the continuos transfer functions====//\n", +"Ts=0.1;//Sampling time ie delta_T\n", +"delay3=3/Ts;\n", +"delay1=1/Ts;\n", +"delay7=7/Ts;\n", +"N=100/Ts;//Model Order\n", +"s=%s;\n", +"G=syslin('c',diag(matrix(G,1,4)));//Transfer function\n", +"t=0:Ts:N*Ts;\n", +"u_sim=ones(4,length(t));\n", +"//Modeling Output delays through input delay in steps\n", +"u_sim(1,1:(delay1))=zeros(1,delay1); \n", +"u_sim(3,1:(delay7))=zeros(1,delay7); \n", +"u_sim([2 4],1:(delay3))=zeros(2,delay3);\n", +"S=csim(u_sim,t,G)';//generating step response model for real plant\n", +"//plot(t,S);\n", +"S(1,:)=[];\n", +"//Now we have these step response models for each of the transfer functions\n", +"//[S1 S3\n", +"//S2 S4 \n", +"T=120;//Simulation Run Time in minutes\n", +"n=T/Ts*2+1; //no. of discrete points in our domain of analysis\n", +"//===============================================================//\n", +"//===============================================================//\n", +"//===============================================================//\n", +"//========Set point as +1 in X-D==X-B loop in manual=============//\n", +"//p is the controller output\n", +"p=zeros(n,2);\n", +"delta_p=zeros(n,2);\n", +"e=zeros(n,2); //errors=(ysp-y) on which PI acts\n", +"ysp=zeros(n,2);\n", +"ysp((n-1)/2+1:n,1)=ones(n-((n-1)/2+1)+1,1);\n", +"t=-(n-1)/2*Ts:Ts:(n-1)/2*Ts;\n", +"y=zeros(n,2);\n", +"for k=(n-1)/2+1:n\n", +" \n", +" //Error e\n", +" e(k,:)=ysp(k-1,:)-y(k-1,:);\n", +" delta_e(k,:)=e(k,:)-e(k-1,:);\n", +" //Controller calculation----Digital PID----Eqn 7-28 Pg 136 (Velocity form)\n", +" p(k,1)=p(k-1,1)+([delta_e(k,1)+e(k,1)*diag(Ts/taui1)]*diag(Kc1));\n", +" //1-1/2-2 pairing\n", +" \n", +" delta_p(k,:)=p(k,:)-p(k-1,:);\n", +" \n", +" //Output\n", +" y(k,1)=[S(1:N-1,1);S(1:N-1,3)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1)]...\n", +" +[S(N,1) S(N,3)]*[p(k-N,1);p(k-N,2)];\n", +" y(k,2)=[S(1:N-1,2);S(1:N-1,4)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);]...\n", +" +[S(N,2) S(N,4)]*[p(k-N,1);p(k-N,2)];\n", +"end\n", +"subplot(2,1,1);\n", +"plot(t',y(:,1),'b-');\n", +"set(gca(),'data_bounds',[0 40 0 1.4]); //putting bounds on display\n", +"l=legend('$x_B\text{ loop in manual}$',position=1);\n", +"xtitle('','Time(min)','$x_D$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;\n", +"//===============================================================//\n", +"//===============================================================//\n", +"//===============================================================//\n", +"//========Set point as +1 in X-B==X-D loop in manual=============//\n", +"//p is the controller output\n", +"p=zeros(n,2);\n", +"delta_p=zeros(n,2);\n", +"e=zeros(n,2); //errors=(ysp-y) on which PI acts\n", +"ysp=zeros(n,2);\n", +"ysp((n-1)/2+1:n,2)=ones(n-((n-1)/2+1)+1,1);\n", +"t=-(n-1)/2*Ts:Ts:(n-1)/2*Ts;\n", +"y=zeros(n,2);\n", +"for k=(n-1)/2+1:n\n", +" \n", +" //Error e\n", +" e(k,:)=ysp(k-1,:)-y(k-1,:);\n", +" delta_e(k,:)=e(k,:)-e(k-1,:);\n", +" //Controller calculation----Digital PID----Eqn 7-28 Pg 136 (Velocity form)\n", +" p(k,2)=p(k-1,2)+([delta_e(k,2)+e(k,2)*diag(Ts/taui2)]*diag(Kc2));\n", +" //1-1/2-2 pairing\n", +" \n", +" delta_p(k,:)=p(k,:)-p(k-1,:);\n", +" \n", +" //Output\n", +" y(k,1)=[S(1:N-1,1);S(1:N-1,3)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1)]...\n", +" +[S(N,1) S(N,3)]*[p(k-N,1);p(k-N,2)];\n", +" y(k,2)=[S(1:N-1,2);S(1:N-1,4)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);]...\n", +" +[S(N,2) S(N,4)]*[p(k-N,1);p(k-N,2)];\n", +"end\n", +"subplot(2,1,2);\n", +"plot(t',y(:,2),'b-');\n", +"set(gca(),'data_bounds',[0 40 0 1.4]); //putting bounds on display\n", +"l=legend('$x_D\text{ loop in manual}$',position=1);\n", +"xtitle('','Time(min)','$x_B$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;\n", +"//===============================================================//\n", +"//===============================================================//\n", +"//===============================================================//\n", +"//========Set point as +1 in X-D==Both loops Automatic===========//\n", +"//p is the controller output\n", +"p=zeros(n,2);\n", +"delta_p=zeros(n,2);\n", +"e=zeros(n,2); //errors=(ysp-y) on which PI acts\n", +"ysp=zeros(n,2);\n", +"ysp((n-1)/2+1:n,1)=ones(n-((n-1)/2+1)+1,1);\n", +"t=-(n-1)/2*Ts:Ts:(n-1)/2*Ts;\n", +"y=zeros(n,2);\n", +"for k=(n-1)/2+1:n\n", +" \n", +" //Error e\n", +" e(k,:)=ysp(k-1,:)-y(k-1,:);\n", +" delta_e(k,:)=e(k,:)-e(k-1,:);\n", +" //Controller calculation----Digital PID----Eqn 7-28 Pg 136 (Velocity form)\n", +"// p(k,:)=p(k-1,:)+flipdim([delta_e(k,:)+e(k,:)*diag(Ts./tauI)]*diag(Kc),2);\n", +" p(k,:)=p(k-1,:)+([delta_e(k,:)+e(k,:)*diag(Ts./tauI)]*diag(Kc));\n", +" //1-1/2-2 pairing\n", +" \n", +" delta_p(k,:)=p(k,:)-p(k-1,:);\n", +" \n", +" //Output\n", +" y(k,1)=[S(1:N-1,1);S(1:N-1,3)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1)]...\n", +" +[S(N,1) S(N,3)]*[p(k-N,1);p(k-N,2)];\n", +" y(k,2)=[S(1:N-1,2);S(1:N-1,4)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);]...\n", +" +[S(N,2) S(N,4)]*[p(k-N,1);p(k-N,2)];\n", +"end\n", +"subplot(2,1,1);\n", +"plot(t',y(:,1),'r--');\n", +"set(gca(),'data_bounds',[0 40 0 1.4]); //putting bounds on display\n", +"l=legend('$x_B\text{ loop in manual}$','$\text{Both loops in automatic}$',position=1);\n", +"//l.font_size=5;\n", +"xtitle('','Time(min)','$x_D$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;\n", +"//===============================================================//\n", +"//===============================================================//\n", +"//===============================================================//\n", +"//========Set point as +1 in X-B==Both loops Automatic===========//\n", +"//p is the controller output\n", +"p=zeros(n,2);\n", +"delta_p=zeros(n,2);\n", +"e=zeros(n,2); //errors=(ysp-y) on which PI acts\n", +"ysp=zeros(n,2);\n", +"ysp((n-1)/2+1:n,2)=ones(n-((n-1)/2+1)+1,1);\n", +"t=-(n-1)/2*Ts:Ts:(n-1)/2*Ts;\n", +"y=zeros(n,2);\n", +"for k=(n-1)/2+1:n\n", +" \n", +" //Error e\n", +" e(k,:)=ysp(k-1,:)-y(k-1,:);\n", +" delta_e(k,:)=e(k,:)-e(k-1,:);\n", +" //Controller calculation----Digital PID----Eqn 7-28 Pg 136 (Velocity form)\n", +"// p(k,:)=p(k-1,:)+flipdim([delta_e(k,:)+e(k,:)*diag(Ts./tauI)]*diag(Kc),2);\n", +" p(k,:)=p(k-1,:)+([delta_e(k,:)+e(k,:)*diag(Ts./tauI)]*diag(Kc));\n", +" //1-1/2-2 pairing\n", +" \n", +" delta_p(k,:)=p(k,:)-p(k-1,:);\n", +" \n", +" //Output\n", +" y(k,1)=[S(1:N-1,1);S(1:N-1,3)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1)]...\n", +" +[S(N,1) S(N,3)]*[p(k-N,1);p(k-N,2)];\n", +" y(k,2)=[S(1:N-1,2);S(1:N-1,4)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);]...\n", +" +[S(N,2) S(N,4)]*[p(k-N,1);p(k-N,2)];\n", +"end\n", +"subplot(2,1,2);\n", +"plot(t',y(:,2),'r--');\n", +"set(gca(),'data_bounds',[0 40 0 1.4]); //putting bounds on display\n", +"l=legend('$x_D\text{ loop in manual}$','$\text{Both loops in automatic}$',position=1);\n", +"xtitle('','Time(min)','$x_B$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;\n", +"//Also refer to Example 22.4 for similar application of algorithm of multiploop PID" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 16.6: Sensitivity_of_steady_state_gain_matrix.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 16.6\n", +"disp('Example 16.6')\n", +"K1=[1 0;10 1];//K with K12=0\n", +"eig1=spec(K1);\n", +"sigma1=spec(K1'*K1);\n", +"CN1=sqrt(max(sigma1)/min(sigma1))\n", +"mprintf('\nEigenvalues of K1 are %f and %f\n and CN is %f',eig1',CN1)\n", +"K2=[1 0.1;10 1];//K with K12=0.1\n", +"eig2=spec(K2);\n", +"sigma2=spec(K2'*K2);\n", +"CN2=sqrt(max(sigma2)/min(sigma2))\n", +"mprintf('\nEigenvalues of K2 are %f and %f\n and CN is %f',eig2',CN2)" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 16.7: Preferred_multiloop_control_strategy.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 16.7\n", +"disp('Example 16.7')\n", +"X=[0.48 0.9 -0.006;0.52 0.95 0.008; 0.90 -0.95 0.020];\n", +"[U,S,V]=svd(X)\n", +"RGA=X.*([inv(X)]') //Eqn 16-36\n", +"//Condition no. of X\n", +"CN=max(diag(S))/min(diag(S))\n", +"//Note that condition no. can also be found with command cond(X)\n", +"// The RGA given in the book is wrong! Eqn 16-73 is wrong.\n", +"mprintf('\n The RGA given in the book is wrong! Eqn 16-73 is wrong.\n')\n", +"disp(RGA,'RGA=')\n", +"X1=X(1:2,1:2);\n", +"X2=X(1:2,[1 3]);\n", +"X3=X(1:2,2:3);\n", +"X4=X([1 3],1:2);\n", +"X5=X([1 3],[1 3]);\n", +"X6=X([1 3],2:3);\n", +"X7=X([2 3],1:2);\n", +"X8=X([2 3],[1 3]);\n", +"X9=X([2 3],2:3);\n", +"lamda1=max(X1.*inv(X1'));\n", +"lamda2=max(X2.*inv(X2'));\n", +"lamda3=max(X3.*inv(X3'));\n", +"lamda4=max(X4.*inv(X4'));\n", +"lamda5=max(X5.*inv(X5'));\n", +"lamda6=max(X6.*inv(X6'));\n", +"lamda7=max(X7.*inv(X7'));\n", +"lamda8=max(X8.*inv(X8'));\n", +"lamda9=max(X9.*inv(X9'));\n", +"mprintf('\n Pairing no. CN lambda\n')\n", +"mprintf('\n 1 %f %f',cond(X1),lamda1)\n", +"mprintf('\n 2 %f %f',cond(X2),lamda2)\n", +"mprintf('\n 3 %f %f',cond(X3),lamda3)\n", +"mprintf('\n 4 %f %f',cond(X4),lamda4)\n", +"mprintf('\n 5 %f %f',cond(X5),lamda5)\n", +"mprintf('\n 6 %f %f',cond(X6),lamda6)\n", +"mprintf('\n 7 %f %f',cond(X7),lamda7)\n", +"mprintf('\n 8 %f %f',cond(X8),lamda8)\n", +"mprintf('\n 9 %f %f',cond(X9),lamda9)" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/17-Digital_Sampling_Filtering_and_Control.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/17-Digital_Sampling_Filtering_and_Control.ipynb new file mode 100644 index 0000000..7d405d5 --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/17-Digital_Sampling_Filtering_and_Control.ipynb @@ -0,0 +1,636 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 17: Digital Sampling Filtering and Control" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 17.1: Performance_of_alternative_filters.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 17.1\n", +"disp('Example 17.1')\n", +"//In this solution we assume that a sampled signal is given to us at a very fast\n", +"//sampling rate and then we resample from it for our computations\n", +"//This depicts how data is in practical situations. \n", +"//Since computers are digital data is always discrete\n", +"//A more kiddish way of writing this code would have been to make a function\n", +"//which takes time as input and gives signal value as output ie generate a \n", +"//continuous signal definition. Then no matter what our sampling time is \n", +"//we can always get the desired values by calling the function say func(Ts*k)\n", +"//where Ts denotes sampling time and k is the index no.(ie deltaT*k)\n", +"//In principle this will also work fine and will reduce the length of the code\n", +"//but this will not lead to learning for coding in practical situations\n", +"Ts=0.001 //sampling time for analog\n", +"t=0:Ts:5;\n", +"n=length(t);\n", +"square_base=0.5*squarewave((t-0.5)*2*%pi/3)+0.5;\n", +"ym=square_base+0.25*sin(t*2*%pi*9);\n", +"subplot(2,2,1)\n", +"plot(t,[square_base' ym'])\n", +"xtitle('Fig 17.6 (a)','Time(min)','Output');\n", +" \n", +"//Analog Filter\n", +"tauf1=0.1;tauf2=0.4;\n", +"s=%s;\n", +"F1=syslin('c',1/(tauf1*s+1));\n", +"F2=syslin('c',1/(tauf2*s+1));\n", +"yf1=csim(ym,t,F1);\n", +"yf2=csim(ym,t,F2);\n", +"subplot(2,2,2);\n", +"plot(t,[yf1' yf2' square_base'])\n", +"legend('$\tau_F=0.1\ min$','$\tau_F=0.4\ min$',position=3);\n", +"xtitle('Fig 17.6 (b)','Time(min)','Output');\n", +" \n", +" //Note that analog filtering can also be achieved by perfect sampling in EWMA digital filter\n", +" //Since Exponentially weighted digital filter is an exact discretization of analog\n", +" //filter if we take Ts=0.001 ie the perfect sampling of data we get identical answers\n", +" //from digital or analog filter. You can try this by chaning Ts1 or Ts2 to 0.001\n", +"//Digital filtering\n", +"Ts1=0.05;Ts2=0.1;\n", +"alpha1=exp(-Ts1/tauf1);\n", +"alpha2=exp(-Ts2/tauf1);\n", +"samples1=1:Ts1/Ts:n;\n", +"samples2=1:Ts2/Ts:n;\n", +"yf1=zeros(length(samples1),1);\n", +"yf2=zeros(length(samples2),1);\n", +"for k=1:length(samples1)-1\n", +" yf1(k+1)=alpha1*yf1(k)+(1-alpha1)*ym(samples1(k));\n", +"end\n", +"for k=1:length(samples2)-1\n", +" yf2(k+1)=alpha2*yf2(k)+(1-alpha2)*ym(samples2(k));\n", +"end\n", +"subplot(2,2,3);\n", +"plot(t(samples1)',[yf1],color='blue');\n", +"plot(t(samples2)',yf2,color='red');\n", +"plot(t,square_base,color='black');\n", +"legend('$\Delta t=0.05 \ min$','$\Delta t=0.1\ min$',position=3);\n", +"xtitle('Fig 17.6 (c)','Time(min)','Output');\n", +" \n", +"//Moving Filter\n", +"N1=3;\n", +"N2=7;\n", +"yf1=zeros(1,length(samples1))\n", +"yf2=zeros(1,length(samples1))\n", +"for k=N1+1:length(samples1)\n", +" yf1(k)=yf1(k-1)+1/N1*(ym(samples1(k))-ym(samples1(k-N1)));\n", +"end\n", +"for k=N2+1:length(samples1)\n", +" yf2(k)=yf2(k-1)+1/N2*(ym(samples1(k))-ym(samples1(k-N2)));\n", +"end\n", +"//for k=N2+1:n\n", +"// yf2(k)=yf2(k-1)+1/N2*(ym(k)-ym(k-N2));\n", +"//end\n", +"subplot(2,2,4);\n", +"plot(t(samples1),[yf1' yf2'])\n", +"plot(t,square_base,color='black');\n", +"legend('$N^*=3$','$N^*=7$',position=4);\n", +"xtitle('Fig 17.6 (d)','Time(min)','Output');\n", +" \n", +"//Now for the gaussian noise\n", +"scf();\n", +"Ts=0.001 //sampling time for analog\n", +"t=0:Ts:5;\n", +"n=length(t);\n", +"square_base=0.5*squarewave((t-0.5)*2*%pi/3)+0.5;\n", +"ym=square_base+grand(1,length(t),'nor', 0, sqrt(0.1));//0.1 is for setting variance=0.1\n", +"subplot(2,2,1)\n", +"plot(t,[square_base' ym'])\n", +"xtitle('Fig 17.6 (a)','Time(min)','Output');\n", +" \n", +"//Analog Filter\n", +"tauf1=0.1;tauf2=0.4;\n", +"s=%s;\n", +"F1=syslin('c',1/(tauf1*s+1));\n", +"F2=syslin('c',1/(tauf2*s+1));\n", +"yf1=csim(ym,t,F1);\n", +"yf2=csim(ym,t,F2);\n", +"subplot(2,2,2);\n", +"plot(t,[yf1' yf2' square_base'])\n", +"legend('$\tau_F=0.1\ min$','$\tau_F=0.4\ min$',position=3);\n", +"xtitle('Fig 17.6 (b)','Time(min)','Output');\n", +" \n", +"//Digital filtering\n", +"Ts1=0.05;Ts2=0.1;\n", +"alpha1=exp(-Ts1/tauf1);\n", +"alpha2=exp(-Ts2/tauf1);\n", +"samples1=1:Ts1/Ts:n;\n", +"samples2=1:Ts2/Ts:n;\n", +"yf1=zeros(length(samples1),1);\n", +"yf2=zeros(length(samples2),1);\n", +"for k=1:length(samples1)-1\n", +" yf1(k+1)=alpha1*yf1(k)+(1-alpha1)*ym(samples1(k));\n", +"end\n", +"for k=1:length(samples2)-1\n", +" yf2(k+1)=alpha2*yf2(k)+(1-alpha2)*ym(samples2(k));\n", +"end\n", +"subplot(2,2,3);\n", +"plot(t(samples1)',[yf1],color='blue');\n", +"plot(t(samples2)',yf2,color='red');\n", +"plot(t,square_base,color='black');\n", +"legend('$\Delta t=0.05 \ min$','$\Delta t=0.1\ min$',position=3);\n", +"xtitle('Fig 17.6 (c)','Time(min)','Output');\n", +" \n", +"//Moving Filter\n", +"N1=3;\n", +"N2=7;\n", +"yf1=zeros(1,length(samples1))\n", +"yf2=zeros(1,length(samples1))\n", +"for k=N1+1:length(samples1)\n", +" yf1(k)=yf1(k-1)+1/N1*(ym(samples1(k))-ym(samples1(k-N1)));\n", +"end\n", +"for k=N2+1:length(samples1)\n", +" yf2(k)=yf2(k-1)+1/N2*(ym(samples1(k))-ym(samples1(k-N2)));\n", +"end\n", +"//for k=N2+1:n\n", +"// yf2(k)=yf2(k-1)+1/N2*(ym(k)-ym(k-N2));\n", +"//end\n", +"subplot(2,2,4);\n", +"plot(t(samples1),[yf1' yf2'])\n", +"plot(t,square_base,color='black');\n", +"legend('$N^*=3$','$N^*=7$',position=4);\n", +"xtitle('Fig 17.6 (d)','Time(min)','Output');\n", +" \n", +"mprintf('Please note that for guassian noise\n results...\n", +" will be different from book owing to randomness\n...\n", +" we do not know the seed for the random noise')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 17.2: Response_of_first_order_difference_equation.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 17.2\n", +"disp('Example 17.2')\n", +"Ts=1;//sampling time\n", +"K=2;\n", +"tau=1;\n", +"alpha=exp(-Ts/tau)\n", +"n=10;\n", +"y=zeros(n,1)\n", +"u=1; //input\n", +"for i=1:n\n", +" y(i+1)=alpha*y(i)+K*(1-alpha)*u;\n", +"end\n", +"disp(y,'yk=')\n", +"mprintf('\n Note that in the book K=20 is wrong, it should be K=2\n...\n", +" that is a first order function with gain 2 is given an input step')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 17.3: Recursive_relation_with_inputs.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 17.3\n", +"disp('Example 17.3')\n", +"z=%z;\n", +"Gz=(-0.3225*z^-2+0.5712*z^-3)/(1-0.9744*z^-1+0.2231*z^-2);\n", +"G=tf2ss(Gz)\n", +"n=10;\n", +"u=ones(1,n);\n", +"y=dsimul(G,u);\n", +"disp(y','y=')\n", +"mprintf('\n\nAlternatively the simulation can also be done\n...\n", +"using syslin(d,Gz) and flts(u,Gz)\n\n')\n", +"Gz2=syslin('d',Gz);\n", +"y2=flts(u,Gz2)\n", +"disp(y2','y2=')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 17.5: Digital_control_of_pressure_in_a_tank.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 17.5\n", +"disp('Example 17.5')\n", +"deltaT=[0.05 0.25 0.5 1]';//sampling time\n", +"K=-20;\n", +"theta=1+(deltaT/2);//Add half of sampling time to delay for finding PI settings\n", +"tau=5;\n", +"//Table 11.3 ITAE disturbance settings\n", +"//Note that there is an error in book solution saying Table 11.2\n", +"//It should be table 11.3\n", +"Y=0.859*(theta/tau)^(-0.977);Kc=Y/K;\n", +"taui=tau*(0.674*(theta/tau)^-0.680).^-1;\n", +"mprintf('\n ITAE(disturbance) \n')\n", +"mprintf(' deltaT Kc tauI')\n", +"mprintf('\n %f %f %f',deltaT,Kc,taui)\n", +"//Finding digital controller settings\n", +"//Eqn 17-55\n", +"a0=1+deltaT./taui;\n", +"a1=-(1); //since tauD=0\n", +"a2=0;\n", +"z=%z;\n", +"Gcz=Kc.*(a0+a1*z^-1)./(1-z^-1);\n", +"//Refer to table 17.1 to convert continuous transfer function to digital form\n", +"Gp=K*(1-exp(-1/tau*deltaT)).*z^(-1+(-1)./deltaT)./(1-exp((-1)/tau*deltaT)*z^-1);//z^(-1/deltaT) for delay\n", +"G_CL=syslin('d',((Gp)./(Gcz.*Gp+1)));\n", +"t=0:deltaT(1):15\n", +"u=ones(1,length(t));\n", +"yt=flts(u,G_CL(1,1));\n", +"plot(t,yt,'-')\n", +"t=0:deltaT(2):15\n", +"u=ones(1,length(t));\n", +"yt=flts(u,G_CL(2,1));\n", +"plot(t,yt,'green--')\n", +"t=0:deltaT(3):15\n", +"u=ones(1,length(t));\n", +"yt=flts(u,G_CL(3,1));\n", +"plot(t,yt,'black-.')\n", +"t=0:deltaT(4):15\n", +"u=ones(1,length(t));\n", +"yt=flts(u,G_CL(4,1));\n", +"plot(t,yt,'red:')\n", +"set(gca(),'data_bounds',[0 15 -8 1]); //putting bounds on display\n", +"l=legend('$\Delta t=0.05\ min$','$\Delta t=0.25\ min$','$\Delta t=0.5\ min$','$\Delta t=1\ min$',position=4);\n", +"xtitle('Example 17.5','Time(min)','$y$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;\n", +"mprintf('\nNote that there is a mismatch between the book simulation and what\n...\n", +"what we get from SCILAB. The book is wrong. This has been crosschecked using\n...\n", +"simulation in SIMULINK(MATLAB)')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 17.6: Dahlin_controller.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 17.6\n", +"disp('Example 17.6')\n", +"//Note that for solving this example there are two ways\n", +"//One is to do this in xcos which is very easy to do\n", +"//and one can learn the same from example 17.5's solution\n", +"//To get the controller outputs at every point in xcos\n", +"//just add a scope to the leg connecting controller and\n", +"//zero order hold unit before the continuous time block\n", +"//The other method is given here so that the reader learns more\n", +"//of what all can be done in scilab\n", +"//Here we deal with the controller in time domain rather than z domain\n", +"z=%z;\n", +"N=0;\n", +"a1=-1.5353;\n", +"a2=0.5866;\n", +"b1=0.0280;\n", +"b2=0.0234;\n", +"G=(b1+b2*z^-1)*z^(-N-1)/(1+a1*z^-1+a2*z^-2);\n", +"h=0;//no process delay\n", +"s=%s;\n", +"lamda=1;\n", +"Y_Ysp=1/(lamda*s+1);//exp(-h*s) is one because h=0 Eqn 17-62\n", +"Ts=1;//sampling time\n", +"A=exp(-Ts/lamda);\n", +"//Eqn 17-63\n", +"Y_Ysp_d=(1-A)*z^(-N-1)/(1-A*z^-1);\n", +"G_DC=1/G*(Y_Ysp_d)/(1-Y_Ysp_d); //Eqn 17-61\n", +"ysp=[zeros(1,4) ones(1,16)]\n", +"Gz_CL=syslin('d',G*G_DC/(G*G_DC+1));//Closed loop discrete system\n", +"yd=flts(ysp,Gz_CL) //Discrete Output due to set point change\n", +"//plot(yd)\n", +"e=ysp-yd; //Since we know set point and the output of the system we can use \n", +"//this info to find out the errors at the discrete time points\n", +"//note that here we have exploited in a very subtle way the property of a \n", +"//discrete system that only the values at discrete points matter for\n", +"//any sort of calculation\n", +"//Now this error can be used to find out the controller effort\n", +"e_coeff=coeff(numer(G_DC));\n", +"p_coeff=coeff(denom(G_DC));\n", +"n=20;//Time in minutes discretized with Ts=1 min\n", +"p=zeros(1,n); //Controller effort\n", +"for k=3:n\n", +" p(k)=(-p_coeff(2)*p(k-1)-p_coeff(1)*p(k-2)+e_coeff*[e(k-2) e(k-1) e(k)]')/p_coeff(3);\n", +"end\n", +"subplot(2,2,2)\n", +"plot2d2(p)\n", +"xtitle('Fig 17.11 (a)','Time(min)','Dahlin Controller effort (p)');\n", +"//Now we simulate the continuous version of the plant to get output in between\n", +"//the discrete point. This will help us ascertain the efficacy of the controller\n", +"//at points other than the discrete points\n", +"//Note that this is required to be checked because deltaT=1. had it been much\n", +"//smaller like 0.01 it would have been a good approx to a continuous system\n", +"//thus making this interpolation check redundant\n", +"s=%s;\n", +"Gp=syslin('c',1/(5*s+1)/(3*s+1));//continuous time version of process \n", +"Ts_c=0.01;//sampling time for continuous system\n", +"t=Ts_c:Ts_c:length([0 p])*Ts;\n", +"p_c=matrix(repmat([0 p],Ts/Ts_c,1),1,Ts/Ts_c*length([0 p]))//hack for zero order hold\n", +"//p_c means controller effort which is continous\n", +"yc=csim(p_c,t,Gp);\n", +"subplot(2,2,1)\n", +"plot(t,yc)\n", +"plot2d2(ysp)\n", +"legend('Dahlin Controller','Set point',position=4)\n", +"xtitle('Fig 17.11 (a)','Time(min)','Output');\n", +"//=============Now we do calculations for modified Dahlin controller========//\n", +"//==========================================================================//\n", +"//Y_Ysp_d=(1-A)*z^(-N-1)/(1-A*z^-1)*(b1+b2*z^-1)/(b1+b2); //Vogel Edgar\n", +"//Page 362 just after solved example\n", +"G_DC_bar=(1-1.5353*z^-1+0.5866*z^-2)/(0.0280+0.0234)*0.632/(1-z^-1); \n", +"//G_DC2=1/G*((1-A)*z^(-N-1))/(1-A*z^-1-(1-A)*z^(-N-1)); //Eqn 17-61\n", +"//G_DC=(1-1.5353*z^-1+0.5866*z^-2)/(0.0280+0.0234*z^-1)*0.632/(1-z^-1);\n", +"ysp=[zeros(1,4) ones(1,16)]\n", +"Gz_CL=syslin('d',G*G_DC_bar/(G*G_DC_bar+1));//Closed loop discrete system\n", +"yd=flts(ysp,Gz_CL) //Discrete Output due to set point change\n", +"//plot(yd)\n", +"e=ysp-yd; //Since we know set point and the output of the system we can use \n", +"//this info to find out the errors at the discrete time points\n", +"//note that here we have exploited in a very subtle way the property of a \n", +"//discrete system that only the values at discrete points matter for\n", +"//any sort of calculation\n", +"//Now this error can be used to find out the controller effort\n", +"e_coeff=coeff(numer(G_DC_bar));\n", +"p_coeff=coeff(denom(G_DC_bar));\n", +"n=20;//Time in minutes discretized with Ts=1 min\n", +"p=zeros(1,n); //Controller effort\n", +"for k=3:n\n", +" p(k)=(-p_coeff(2)*p(k-1)-p_coeff(1)*p(k-2)+e_coeff*[e(k-2) e(k-1) e(k)]')/p_coeff(3);\n", +"end\n", +"subplot(2,2,4)\n", +"plot2d2(p)\n", +"xtitle('Fig 17.11 (b)','Time(min)','Modified Dahlin Controller effort (p)');\n", +"//Now we simulate the continuous version of the plant to get output in between\n", +"//the discrete point. This will help us ascertain the efficacy of the controller\n", +"//at points other than the discrete points\n", +"//Note that this is required to be checked because deltaT=1. had it been much\n", +"//smaller like 0.01 it would have been a good approx to a continuous system\n", +"//thus making this interpolation check redundant\n", +"s=%s;\n", +"Gp=syslin('c',1/(5*s+1)/(3*s+1));//continuous time version of process \n", +"Ts_c=0.01;//sampling time for continuous system\n", +"t=Ts_c:Ts_c:length([0 p])*Ts;\n", +"p_c=matrix(repmat([0 p],Ts/Ts_c,1),1,Ts/Ts_c*length([0 p]))//hack for zero order hold\n", +"//p_c means controller effort which is continous\n", +"yc=csim(p_c,t,Gp);\n", +"subplot(2,2,3)\n", +"plot(t,yc)\n", +"plot2d2(ysp)\n", +"legend('Modified Dahlin Controller','Set point',position=4)\n", +"xtitle('Fig 17.11 (b)','Time(min)','Output');" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 17.7: Non_ringing_Dahlin_controller.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 17.7\n", +"disp('Example 17.7')\n", +"//Note that for solving this example there are two ways\n", +"//One is to do this in xcos which is very easy to do\n", +"//and one can learn the same from example 17.5's solution\n", +"//To get the controller outputs at every point in xcos\n", +"//just add a scope to the leg connecting controller and\n", +"//zero order hold unit before the continuous time block\n", +"//The other method is given here so that the reader learns more\n", +"//of what all can be done in scilab\n", +"//Here we deal with the controller in time domain rather than z domain\n", +"z=%z;\n", +"N=0;\n", +"a1=-1.5353;\n", +"a2=0.5866;\n", +"b1=0.0280;\n", +"b2=0.0234;\n", +"G=(b1+b2*z^-1)*z^(-N-1)/(1+a1*z^-1+a2*z^-2);\n", +"h=0;//no process delay\n", +"s=%s;\n", +"lamda=1;\n", +"Y_Ysp=1/(lamda*s+1);//exp(-h*s) is one because h=0 Eqn 17-62\n", +"Ts=1;//sampling time\n", +"A=exp(-Ts/lamda);\n", +"//=============Now we do calculations for modified Dahlin controller========//\n", +"//==========================================================================//\n", +"//Page 362 just after solved example\n", +"G_DC_bar=(1-1.5353*z^-1+0.5866*z^-2)/(0.0280+0.0234)*0.632/(1-z^-1); \n", +"ysp=[zeros(1,4) ones(1,16)]\n", +"Gz_CL=syslin('d',G*G_DC_bar/(G*G_DC_bar+1));//Closed loop discrete system\n", +"yd=flts(ysp,Gz_CL) //Discrete Output due to set point change\n", +"//plot(yd)\n", +"e=ysp-yd; //Since we know set point and the output of the system we can use \n", +"//this info to find out the errors at the discrete time points\n", +"//note that here we have exploited in a very subtle way the property of a \n", +"//discrete system that only the values at discrete points matter for\n", +"//any sort of calculation\n", +"//Now this error can be used to find out the controller effort\n", +"e_coeff=coeff(numer(G_DC_bar));\n", +"p_coeff=coeff(denom(G_DC_bar));\n", +"n=20;//Time in minutes discretized with Ts=1 min\n", +"p=zeros(1,n); //Controller effort\n", +"for k=3:n\n", +" p(k)=(-p_coeff(2)*p(k-1)-p_coeff(1)*p(k-2)+e_coeff*[e(k-2) e(k-1) e(k)]')/p_coeff(3);\n", +"end\n", +"subplot(3,2,2)\n", +"plot2d2(p)\n", +"xtitle('Fig 17.12 (a)','Time(min)','Modified Dahlin Controller effort (p)');\n", +"//Now we simulate the continuous version of the plant to get output in between\n", +"//the discrete point. This will help us ascertain the efficacy of the controller\n", +"//at points other than the discrete points\n", +"//Note that this is required to be checked because deltaT=1. had it been much\n", +"//smaller like 0.01 it would have been a good approx to a continuous system\n", +"//thus making this interpolation check redundant\n", +"s=%s;\n", +"Gp=syslin('c',1/(5*s+1)/(3*s+1));//continuous time version of process \n", +"Ts_c=0.01;//sampling time for continuous system\n", +"t=Ts_c:Ts_c:length([0 p])*Ts;\n", +"p_c=matrix(repmat([0 p],Ts/Ts_c,1),1,Ts/Ts_c*length([0 p]))//hack for zero order hold\n", +"//p_c means controller effort which is continous\n", +"yc=csim(p_c,t,Gp);\n", +"subplot(3,2,1)\n", +"plot(t,yc)\n", +"plot2d2(ysp)\n", +"legend('Modified Dahlin Controller','Set point',position=4)\n", +"xtitle('Fig 17.12 (a)','Time(min)','Output');\n", +"//=============Now we do calculations for PID-BD controller========//\n", +"//==========================================================================//\n", +"G_BD=4.1111*(3.1486-5.0541*z^-1+2.0270*z^-2)/(1.7272-2.4444*z^-1+0.7222*z^-2)\n", +"ysp=[zeros(1,4) ones(1,16)]\n", +"Gz_CL=syslin('d',G*G_BD/(G*G_BD+1));//Closed loop discrete system\n", +"yd=flts(ysp,Gz_CL) //Discrete Output due to set point change\n", +"//plot(yd)\n", +"e=ysp-yd; //Since we know set point and the output of the system we can use \n", +"//this info to find out the errors at the discrete time points\n", +"//note that here we have exploited in a very subtle way the property of a \n", +"//discrete system that only the values at discrete points matter for\n", +"//any sort of calculation\n", +"//Now this error can be used to find out the controller effort\n", +"e_coeff=coeff(numer(G_BD));\n", +"p_coeff=coeff(denom(G_BD));\n", +"n=20;//Time in minutes discretized with Ts=1 min\n", +"p=zeros(1,n); //Controller effort\n", +"for k=3:n\n", +" p(k)=(-p_coeff(2)*p(k-1)-p_coeff(1)*p(k-2)+e_coeff*[e(k-2) e(k-1) e(k)]')/p_coeff(3);\n", +"end\n", +"subplot(3,2,4)\n", +"plot2d2(p)\n", +"xtitle('Fig 17.12 (b)','Time(min)','BD Controller effort (p)');\n", +"//Now we simulate the continuous version of the plant to get output in between\n", +"//the discrete point. This will help us ascertain the efficacy of the controller\n", +"//at points other than the discrete points\n", +"//Note that this is required to be checked because deltaT=1. had it been much\n", +"//smaller like 0.01 it would have been a good approx to a continuous system\n", +"//thus making this interpolation check redundant\n", +"s=%s;\n", +"Gp=syslin('c',1/(5*s+1)/(3*s+1));//continuous time version of process \n", +"Ts_c=0.01;//sampling time for continuous system\n", +"t=Ts_c:Ts_c:length([0 p])*Ts;\n", +"p_c=matrix(repmat([0 p],Ts/Ts_c,1),1,Ts/Ts_c*length([0 p]))//hack for zero order hold\n", +"//p_c means controller effort which is continous\n", +"yc=csim(p_c,t,Gp);\n", +"subplot(3,2,3)\n", +"plot(t,yc)\n", +"plot2d2(ysp)\n", +"legend('PID-BD Controller','Set point',position=4)\n", +"xtitle('Fig 17.12 (b)','Time(min)','Output');\n", +"//=============Now we do calculations for Vogel Edgar Dahlin controller========//\n", +"//==========================================================================//\n", +"Y_Ysp_d=(1-A)*z^(-N-1)/(1-A*z^-1)*(b1+b2*z^-1)/(b1+b2); //Vogel Edgar Eqn 17-70\n", +"G_VE=1/G*(Y_Ysp_d)/(1-Y_Ysp_d); //Eqn 17-61\n", +"ysp=[zeros(1,4) ones(1,16)]\n", +"Gz_CL=syslin('d',G*G_VE/(G*G_VE+1));//Closed loop discrete system\n", +"yd=flts(ysp,Gz_CL) //Discrete Output due to set point change\n", +"//plot(yd)\n", +"e=ysp-yd; //Since we know set point and the output of the system we can use \n", +"//this info to find out the errors at the discrete time points\n", +"//note that here we have exploited in a very subtle way the property of a \n", +"//discrete system that only the values at discrete points matter for\n", +"//any sort of calculation\n", +"//Now this error can be used to find out the controller effort\n", +"e_coeff=coeff(numer(G_VE));\n", +"p_coeff=coeff(denom(G_VE));\n", +"n=20;//Time in minutes discretized with Ts=1 min\n", +"p=zeros(1,n); //Controller effort\n", +"for k=3:n\n", +" p(k)=(-p_coeff(2)*p(k-1)-p_coeff(1)*p(k-2)+e_coeff*[e(k-2) e(k-1) e(k)]')/p_coeff(3);\n", +"end\n", +"subplot(3,2,6)\n", +"plot2d2(p)\n", +"xtitle('Fig 17.12 (c)','Time(min)','Vogel Edgar Controller effort (p)');\n", +"//Now we simulate the continuous version of the plant to get output in between\n", +"//the discrete point. This will help us ascertain the efficacy of the controller\n", +"//at points other than the discrete points\n", +"//Note that this is required to be checked because deltaT=1. had it been much\n", +"//smaller like 0.01 it would have been a good approx to a continuous system\n", +"//thus making this interpolation check redundant\n", +"s=%s;\n", +"Gp=syslin('c',1/(5*s+1)/(3*s+1));//continuous time version of process \n", +"Ts_c=0.01;//sampling time for continuous system\n", +"t=Ts_c:Ts_c:length([0 p])*Ts;\n", +"p_c=matrix(repmat([0 p],Ts/Ts_c,1),1,Ts/Ts_c*length([0 p]))//hack for zero order hold\n", +"//p_c means controller effort which is continous\n", +"yc=csim(p_c,t,Gp);\n", +"subplot(3,2,5)\n", +"plot(t,yc)\n", +"plot2d2(ysp)\n", +"legend('Vogel Edgar Controller','Set point',position=4)\n", +"xtitle('Fig 17.12 (c)','Time(min)','Output');\n", +"mprintf('Note that there is some very slight difference between the \n...\n", +"curves shown in book and that obtained from scilab\n...\n", +"this is simply because of more detailed calculation in scilab ')" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/19-Real_Time_Optimization.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/19-Real_Time_Optimization.ipynb new file mode 100644 index 0000000..d96b230 --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/19-Real_Time_Optimization.ipynb @@ -0,0 +1,169 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 19: Real Time Optimization" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 19.2: Nitration_of_Decane.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 19.2\n", +"disp('Example 19.2')\n", +"function y=f_DNO3(r1)\n", +" D1=0.5;D2=0.5;\n", +" r2=0.4-0.5*r1;\n", +" y=r1*D1/(1+r1)^2/(1+r2)+r2*D2/(1+r1)/(1+r2)^2\n", +"endfunction\n", +"function [f, g, ind] = costf(x, ind)\n", +" f=-f_DNO3(x);//cost is negative of function to be maximised\n", +" g=-derivative(f_DNO3,x);//derivative of the cost function\n", +"endfunction\n", +"[fopt, xopt] = optim(costf,0.5);\n", +"disp(xopt,'Optimum value of r1=')\n", +"disp(-fopt,'Max value of DNO3=')\n", +"mprintf('Note that the answer in book is not as accurate as the one\n...\n", +"calculated from scilab')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 19.3: Refinery_blending_and_production.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 19.3\n", +"disp('Example 19.3')\n", +"//function for minimization\n", +"c=-[-24.5 -16 36 24 21 10]';\n", +"//Equality Constraints\n", +"Aeq=[0.80 0.44 -1 0 0 0;0.05 0.1 0 -1 0 0;0.1 0.36 0 0 -1 0;0.05 0.1 0 0 0 -1];\n", +"beq=zeros(4,1);\n", +"//Inequality Constraints\n", +"A=[0 0 1 0 0 0;0 0 0 1 0 0;0 0 0 0 1 0];\n", +"b=[24000 2000 6000]';\n", +"//Lower bound on x\n", +"lb=zeros(6,1);\n", +"//Initial guess: such that it satisfies Aeq*x0=beq\n", +"x0=zeros(6,1);\n", +"x0(1:2)=[5000 3000]';//Initial guess for x1 and x2\n", +"x0(3:6)=Aeq(:,3:6)\(beq-Aeq(:,1:2)*x0(1:2));//solution of linear equations\n", +"//Note that x0 should also satisfy A*x0<b and lb\n", +"[xopt,fopt]=karmarkar(Aeq,beq,c,x0,[],[],[],[],A,b,lb)\n", +"disp(xopt,'Optimum value of x=')\n", +"mprintf('\nMax value of f=$ %f /day\n',-fopt)\n", +"mprintf('\n Note that the answer in book is not as accurate as the one\n...\n", +" calculated from scilab')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 19.4: Fuel_cost_in_boiler_house.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 19.4\n", +"disp('Example 19.4')\n", +"//Here we have Nonlinear programming problem hence we use optim function\n", +"//Since optim does not have the ability to handle constraints\n", +"//we use the penalty method for optimization\n", +"//ie we make the constraints a part of the cost function such that \n", +"//cost function increases severly for any violation of constraints\n", +"//MATLAB users must be familiar with fmincon function in MATLAB\n", +"//Unfortunately a similar function in Scilab is not yet available\n", +"//Fmincon toolbox development for scilab is under development/testing\n", +"x0=[2 4 4 1]'; //Initial guess\n", +"function y=func(x) //x is 4*1 vector\n", +" P1=4.5*x(1)+0.1*x(1)^2+4*x(2)+0.06*x(2)^2;\n", +" P2=4*x(3)+0.05*x(3)^2+3.5*x(4)+0.2*x(4)^2;\n", +" if (P1>30) then\n", +" c1=abs(P1-30)^2;\n", +" elseif P1<18\n", +" c1=abs(P1-18)^2;\n", +" else c1=0;\n", +" end\n", +" if (P2>25) then\n", +" c2=abs(P2-30)^2;\n", +" elseif P2<14\n", +" c2=abs(P2-18)^2;\n", +" else c2=0;\n", +" end\n", +" c3=abs(P1+P2-50)^2;\n", +" c4=abs(x(2)+x(4)-5)^2;\n", +" y=(x(1)+x(3))+100*(c1+c2+c3+c4);\n", +"endfunction\n", +"function [f, g, ind] = costf(x, ind)\n", +" f=func(x);//cost is negative of function to be maximised\n", +" g=derivative(func,x);//derivative of the cost function\n", +"endfunction\n", +"[fopt, xopt] = optim(costf,'b',zeros(4,1), 10*ones(4,1),x0);\n", +"// 'b',zeros(4,1), 10*ones(4,1) stands for lower and upper bounds on x\n", +"disp(xopt,'Optimum value of x=')\n", +"disp(fopt,'Min value of f=')\n", +"mprintf('Note that the answer in book is not as accurate as the one\n...\n", +"calculated from scilab')" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/2-Theoretical_Models_of_Chemical_Processes.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/2-Theoretical_Models_of_Chemical_Processes.ipynb new file mode 100644 index 0000000..2a6fd3e --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/2-Theoretical_Models_of_Chemical_Processes.ipynb @@ -0,0 +1,252 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 2: Theoretical Models of Chemical Processes" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 2.1: Stirred_tank_blending_process.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 2.1\n", +"disp('Example 2.1')\n", +"w1bar=500;\n", +"w2bar=200;\n", +"x1bar=0.4;\n", +"x2bar=0.75;\n", +"wbar=w1bar+w2bar;\n", +"t=0:0.1:25; //Time scale for plotting of graphs\n", +"//(a)\n", +"xbar=(w1bar*x1bar+w2bar*x2bar)/wbar;\n", +"printf('\n (a) The steady state concentration is %f \n',xbar)\n", +"//(b)\n", +"w1bar=400; //flow rate changes, rest remains same\n", +"wbar=w1bar+w2bar;\n", +"tau=3;\n", +"x0=0.5; \n", +"Cstarb=(w1bar*x1bar+w2bar*x2bar)/wbar; //C* variable\n", +"printf('\n (b) The value of C* is %f',Cstarb)\n", +"printf('\n x(t)=0.5exp(-t/3)+%f(1-exp(-t/3)) \n',Cstarb);\n", +"xtd=0.5*exp(-t/3)+Cstarb*(1-exp(-t/3));\n", +"xtb=0.5*exp(-t/3)+Cstarb*(1-exp(-t/3)); //x(t) for part (b)\n", +"//(c)\n", +"w1bar=500;w2bar=100; //flow rate changes, rest remains same\n", +"wbar=w1bar+w2bar;\n", +"tau=3;\n", +"x0=0.5; \n", +"Cstarc=(w1bar*x1bar+w2bar*x2bar)/wbar; //C* variable\n", +"printf('\n (c) The value of C* is %f',Cstarc)\n", +"printf('\n x(t)=0.5exp(-t/3)+%f(1-exp(-t/3)) \n',Cstarc);\n", +"xtc=0.5*exp(-t/3)+Cstarc*(1-exp(-t/3));\n", +"//(d)\n", +"w1bar=500;w2bar=100;x1bar=0.6;x2bar=0.75; //flow rate changes, rest remains same\n", +"wbar=w1bar+w2bar;\n", +"tau=3;\n", +"x0=0.5; \n", +"Cstard=(w1bar*x1bar+w2bar*x2bar)/wbar; //C* variable\n", +"printf('\n (d) The value of C* is %f',Cstard)\n", +"printf('\n x(t)=0.5exp(-t/3)+%f(1-exp(-t/3)) \n',Cstard);\n", +"xtd=0.5*exp(-t/3)+Cstard*(1-exp(-t/3));\n", +"plot2d(t,[xtd',xtb',xtc'])\n", +"xtitle('Parts b through d','Time(min)','$x(t)$');\n", +"a=legend('$(d)$','$(b)$','$(c)$',position=1);\n", +"a.font_size=5;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;\n", +"//(e)\n", +"xNb=(xtb-x0)/(Cstarb-x0); //Normalized response for part b\n", +"xNc=(xtc-x0)/(Cstarc-x0); //Normalized response for part c\n", +"xNd=(xtd-x0)/(Cstard-x0); //Normalized response for part d\n", +"scf() //Creates new window for plotting\n", +"plot2d(t,[xNd',xNb',xNc'],style=[1 1 1]) \n", +"//Style sets the color, -ve values means discrete plotting, +ve means color\n", +"xtitle('Part e','Time(min)','Normalized response');\n", +"a=legend('$(e)$',position=1);\n", +"a.font_size=5;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 2.2: Degrees_of_freedom_1.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 2.2\n", +"disp('Example 2.2')\n", +"N_V=4;\n", +"N_E=1;\n", +"N_F=N_V-N_E;\n", +"printf('\n Degrees of freedom N_F= %i \n',N_F)" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 2.3: Degrees_of_freedom_2.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 2.3\n", +"disp('Example 2.3')\n", +"N_V=7;\n", +"N_E=2;\n", +"N_F=N_V-N_E;\n", +"printf('\n Degrees of freedom N_F= %i \n',N_F)" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 2.4: Electrically_heated_stirred_tank_process.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 2.4\n", +"disp('Example 2.4')\n", +"mprintf('\n Important Note: Errata for book: Values of the parameters \n...\n", +" meCe/heAe and meCe/wC should be 1 min each and not 0.5 min %s \n','')\n", +"Tibar=100;//deg C\n", +"Qbar=5000;//kcal/min\n", +"wc_inv=0.05;// 1/wc degC min/kcal\n", +"//(a)\n", +"Tbar=Tibar+wc_inv*Qbar;\n", +"mprintf('\n (a) Nominal steady state temperature= %i',Tbar)\n", +"mprintf(' degree celsius %s \n','')\n", +"//(b)\n", +"mprintf('\n Eqn 2-29 becomes 10 d2T/dt2 + 12 dT/dt + T = 370 with T(0)=350 %s \n','')\n", +"t=0:0.1:80; //Time values\n", +"Tt_2=350+20*(1-1.089*exp(-t/11.099)+0.084*exp(-t/0.901));//T(t) from order 2 equation\n", +"//(c)\n", +"mprintf('\n Eqn 2-29 becomes 12 dT/dt + T = 370 with T(0)=350 %s \n','')\n", +"Tt_1=350+20*(1-exp(-t/12));//T(t) from order 1 equation\n", +"plot2d(t,[Tt_2',Tt_1'],[2 5],rect=[0 350 80 370])\n", +"xtitle('Ex-2.4','Time(min)','$T(^0C)$');\n", +"a=legend('a Second order','b First order',position=4);\n", +"a.font_size=5;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 2.5: Nonlinear_dynamic_behavior_of_CSTR.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 2.5\n", +"disp('Example 2.5')\n", +"function ydot=CSTR(t,y,Tc) //y is [Conc Temp]' Tc is coolant temp\n", +" q=100;ci=1;V=100;rho=1000;C=0.239;deltaHR=5E4;k0=7.2E10;UA=5E4;Er=8750;\n", +" Ti=350;\n", +" c=y(1);T=y(2);\n", +" k=k0*exp(-Er/T);//Er=E/R\n", +" ydot(1)=1/V*(q*(ci-c)-V*k*c); //ydot(1) is dc_dt\n", +" ydot(2)=1/(V*rho*C)*(q*rho*C*(Ti-T)+deltaHR*V*k*c+UA*(Tc-T))//ydot(2) is dT_dt\n", +"endfunction\n", +"c0=0.5;T0=350;\n", +"y0=[c0 T0]';\n", +"t0=0;\n", +"t=0:0.01:10;\n", +"Tc=[290 305];\n", +"y1 = ode(y0,t0,t,list(CSTR,Tc(1)));\n", +"y2 = ode(y0,t0,t,list(CSTR,Tc(2)));\n", +"y3=[0.5 0;0 350]*ones(2,length(t))\n", +"//Temp plot\n", +"subplot(2,1,1);\n", +"plot(t,[y1(2,:)' y3(2,:)' y2(2,:)']);\n", +"xtitle('Fig 2.7','Time(min)','Reactor Temp(K)');\n", +"legend('290 K','300 K','305 K')\n", +"//conc plot\n", +"subplot(2,1,2);\n", +"plot(t,[y1(1,:)' y3(1,:)' y2(1,:)']);\n", +"xtitle('Fig 2.8','Time(min)','Reactor conc(mol/L)');\n", +"legend('290 K','300 K','305 K');" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/20-Model_Predictive_Control.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/20-Model_Predictive_Control.ipynb new file mode 100644 index 0000000..343ec86 --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/20-Model_Predictive_Control.ipynb @@ -0,0 +1,596 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 20: Model Predictive Control" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 20.1: Step_response_coefficients.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 20.1\n", +"disp('Example 20.1')\n", +"K=5;\n", +"tau=15;//min\n", +"theta=2;//min\n", +"Ts=1;//Sampling period\n", +"k=[0:79]';//samples\n", +"N=80;\n", +"//From eqn 20-5\n", +"S=zeros(N,1);\n", +"S=K*(1-exp(-(k*Ts-theta)/tau));\n", +"S(1:(theta+1),1)=0;//delay\n", +"//Step change\n", +"M=3;\n", +"//Calculating step response from Eqn 20-4\n", +"step=3;//step change occurs at t=3 min\n", +"i=[(theta+1):90]';\n", +"yi=[zeros(theta+step,1);K*M*(1-exp(-(i*Ts-theta)/tau));]\n", +"plot2d(yi,style=-4);\n", +"xtitle('$Ex\ 20.1$','Time(min)','y')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 20.3: J_step_ahead_predictio.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 20.3\n", +"disp('Example 20.3')\n", +"for J=[3 4 6 8] //Tuning parameter\n", +" \n", +"Ts=5;//Sampling time ie delta_T\n", +"N=16;//Model Order\n", +"s=%s;\n", +"G=syslin('c',1/(5*s+1)^5);//Transfer function\n", +"t=0:Ts:N*Ts;\n", +"S=csim('step',t,G)';//generating step response model\n", +"//plot(t,S);\n", +"S(1)=[];\n", +"T=80;//simulation time\n", +"n=T/Ts*2+1;\n", +"u=zeros(1,n); \n", +"//Input initialization 80 min is the Time for simulation\n", +"//We take a few u values in negative time to facilitate usage of step response\n", +"//model\n", +"delta_u=[0 diff(u)];\n", +"yhat_u=zeros(n,1);\n", +"ysp=1;\n", +"for k=(n-1)/2+1+1:n-J //an additional +1 is because MPC starts after set point change\n", +" yhat_u(k+J)=delta_u(k+J-N+1:k-1)*flipdim(S(J+1:N-1),1)+S(N)*u(k+J-N);//unforced predicted y\n", +" disp(yhat_u(k+J))\n", +" delta_u(k)=(ysp-yhat_u(k+J))/S(J);\n", +" u(k)=u(k-1)+delta_u(k);\n", +"end\n", +"u(n-J+1:$)=u(k)*ones(1,J);//Carry forward the u as constant \n", +"t=-(n-1)/2*Ts:Ts:(n-1)*Ts/2;\n", +"subplot(2,1,2);\n", +"if J==3 | J==6 then\n", +" plot2d2(t((n-1)/2+1:n),u((n-1)/2+1:n),style=J);\n", +"end\n", +"legend('J=3','J=6',position=4)\n", +"xtitle('','Time(min)','$u$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;\n", +"res=Ts;//resolution\n", +"//u_cont=matrix(repmat([0 u],res,1),1,res*length([0 u]));\n", +"u_cont=matrix(repmat([u],res,1),1,res*length([u]));\n", +"entries=length(u_cont);\n", +"t_cont=linspace(-T,T+Ts-1,entries);\n", +"yt=csim(u_cont,t_cont,G);\n", +"subplot(2,1,1);\n", +"if J=8 then //for color of plot2d\n", +" J=9\n", +"end\n", +"plot2d(t_cont((entries-Ts)/2+1:$),yt((entries-Ts)/2+1:$),style=J,rect=[0,0,80,1.2]);\n", +"end\n", +"//Other niceties for plots\n", +"subplot(2,1,1);\n", +"plot(t_cont((entries-Ts)/2+1:$),ones(length(t_cont((entries-Ts)/2+1:$)),1),'--');\n", +"legend('J=3','J=4','J=6','J=8',position=4)\n", +"xtitle('Example 20.3','Time(min)','$y$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 20.4: Output_feedback_and_bias_correction.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 20.4\n", +"disp('Example 20.4')\n", +"J=15;\n", +"Ts=1;//Sampling time ie delta_T\n", +"N=80;//Model Order\n", +"s=%s;\n", +"G=syslin('c',5/(15*s+1));//Transfer function\n", +"t=0:Ts:N*Ts;\n", +"u_sim=ones(1,length(t));\n", +"u_sim(1:3)=[0 0 0]; //input delay to account for 2 min delay in G\n", +"S=csim(u_sim,t,G)';//generating step response model for real plant\n", +"//plot(t,S);\n", +"S(1)=[];\n", +"T=100;//simulation time\n", +"n=T/Ts*2+1; //no. of discrete points in our domain of analysis\n", +"//Input initialization T min is the Time for simulation\n", +"//We take a few u values in negative time to facilitate \n", +"//usage of step response model\n", +"u=zeros(n,1);\n", +"d=zeros(n,1);\n", +"delta_u=zeros(n,1);\n", +"delta_u(101+2)=1; //Step change at t=2 min\n", +"u=cumsum(delta_u);\n", +"delta_d=zeros(n,1);\n", +"delta_d(101+8)=0.15;//disturbance t=8 min\n", +"d=cumsum(delta_d);\n", +"yhat=zeros(n,1); //J step ahead predictions\n", +"ytilda=zeros(n,1); //J step ahead predictions corrected\n", +"b=zeros(n,1); //corrections\n", +"t=-(n-1)/2:Ts:(n-1)/2;\n", +"for k=(n-1)/2+1-J:n-J\n", +" yhat(k+J)=S(1:N-1)'*flipdim(delta_u(k+J-N+1:k+J-1),1)+S(N)*u(k+J-N);\n", +" //Predicted y Eqn 20-10\n", +" y(k+J)=S(1:N-1)'*flipdim(delta_u(k+J-N+1:k+J-1),1)+S(N)*u(k+J-N)+...\n", +" S(1:N-1)'*flipdim(delta_d(k+J-N+1:k+J-1),1)+S(N)*d(k+J-N);\n", +" //Actual values of the process\n", +" b(k+J)=y(k)-yhat(k); //Note that there is a delay in corrections\n", +" //which is opposite of prediction\n", +"end\n", +"ytilda=b+yhat;//calculation of corrected y\n", +"plot(t,y,'-',t,yhat,'-.',t,ytilda,'--');\n", +"set(gca(),'data_bounds',[0 100 0 6]); //putting bounds on display\n", +"l=legend('y','$\hat y$','$\tilde y$',position=4);\n", +"l.font_size=5;\n", +"xtitle('Example 20.4','Time(min)','$y$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;\n", +"//======Part b=====//\n", +"G2=syslin('c',4/(20*s+1));//Transfer function\n", +"t2=0:Ts:N*Ts;\n", +"u_sim=ones(1,length(t2));\n", +"u_sim(1:3)=[0 0 0]; //input delay to account for 2 min delay in G\n", +"S2=csim(u_sim,t2,G2)';//generating step response model for model\n", +"//plot(t2,S);\n", +"S2(1)=[];\n", +"yhat=zeros(n,1); //J step ahead predictions\n", +"ytilda=zeros(n,1); //J step ahead predictions corrected\n", +"b=zeros(n,1); //corrections\n", +"for k=(n-1)/2+1-J:n-J\n", +" yhat(k+J)=S2(1:N-1)'*flipdim(delta_u(k+J-N+1:k+J-1),1)+S2(N)*u(k+J-N);\n", +" //Predicted y Eqn 20-10\n", +" y(k+J)=S(1:N-1)'*flipdim(delta_u(k+J-N+1:k+J-1),1)+S(N)*u(k+J-N)+...\n", +" S(1:N-1)'*flipdim(delta_d(k+J-N+1:k+J-1),1)+S(N)*d(k+J-N);\n", +" //Actual values of the process\n", +" b(k+J)=y(k)-yhat(k); //Note that there is a delay in corrections\n", +" //which is opposite of prediction\n", +"end\n", +"ytilda=b+yhat;//calculation of corrected y\n", +"scf();\n", +"plot(t,y,'-',t,yhat,'-.',t,ytilda,'--');\n", +"set(gca(),'data_bounds',[0 100 0 6]); //putting bounds on display\n", +"l=legend('y','$\hat y$','$\tilde y$',position=4);\n", +"l.font_size=5;\n", +"xtitle('Example 20.4','Time(min)','$y$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 20.5: Comparison_of_MPCs_and_PID.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 20.5\n", +"disp('Example 20.5')\n", +"//=============Part(a)========//\n", +"//=============Part(a)========//\n", +"//=============Part(a)========//\n", +"//=============Part(a)========//\n", +"//=============Part(a)========//\n", +"Ts=1;//Sampling time ie delta_T\n", +"N=70;//Model Order\n", +"s=%s;\n", +"G=syslin('c',1/(5*s+1)/(10*s+1));//Transfer function\n", +"t=0:Ts:(N-1)*Ts;\n", +"u_sim=ones(1,length(t));\n", +"//There is automatically an input delay of one unit in csim function\n", +"S=csim(u_sim,t,G)';//generating step response model for real plant\n", +"//plot(t,S);\n", +"T=80;//simulation time\n", +"//Let the three simulations correspond to\n", +"//MPC1==> P=3,M=1\n", +"//MPC2==> P=4,M=2\n", +"///PID==> The PID controller\n", +"P1=3;M1=1;\n", +"P2=4;M2=2;\n", +"S1=S(1:P1);//MPC-1\n", +"S2=[S(1:P2) [0;S(1:P2-1)]];//MPC-2\n", +"//SISO system\n", +"Q=1;\n", +"R=0;//No move suppression\n", +"Kc1=inv(S1'*Q*S1+R*eye(M1,M1))*S1'*Q;//Eqn20-57 MPC1\n", +"Kc2=inv(S2'*Q*S2+R*eye(M2,M2))*S2'*Q;//Eqn20-57 MPC2\n", +"mprintf('\nFor P=3 and M=1, \nKc=\n [%f %f %f]\n',Kc1)\n", +"mprintf('\nFor P=4 and M=2,\nKc=')\n", +"disp(Kc2)\n", +"//=============Part(b)========//\n", +"//=============Part(b)========//\n", +"//=============Part(b)========//\n", +"//=============Part(b)========//\n", +"//=============Part(b)========//\n", +"//=============Part(b)========//\n", +"//=============Part(b) MPC-1========//\n", +"n=T/Ts*2+1; //no. of discrete points in our domain of analysis\n", +"//Input initialization T is the Time for simulation\n", +"//We take a few u values in negative time to facilitate \n", +"//usage of step response model\n", +"u=zeros(n,1);\n", +"d=zeros(n,1);\n", +"delta_u=zeros(n,1);\n", +"//delta_u(101+2)=1; //Step change at t=2 min\n", +"u=cumsum(delta_u);\n", +"delta_d=zeros(n,1);\n", +"//delta_d(101+8)=0.15;//disturbance t=8 min\n", +"d=cumsum(delta_d);\n", +"y=zeros(1,n);//Actual values\n", +"yhat=zeros(1,n); //predicted value\n", +"ydot=zeros(P1,n); //Unforced predictions\n", +"ydottilde=zeros(P1,n); //Corrected unforced predictions\n", +"yr=ones(P1,n);//reference trajectory(same as setpoint)\n", +"edot=zeros(P1,n);//predicted unforced error\n", +"t=-(n-1)/2:Ts:(n-1)/2;\n", +"for k=(n-1)/2+1:n-P1\n", +" \n", +" //Unforced predictions\n", +" for J=1:P1\n", +" ydot(J,k+1)=S(J+1:N-1)'*flipdim(delta_u(k+J-N+1:k-1),1)+S(N)*u(k+J-N);\n", +" end\n", +" //Actual values of the process\n", +" J=0;\n", +" y(k+J)=S(1:N-1)'*flipdim(delta_u(k+J-N+1:k+J-1),1)+S(N)*u(k+J-N)+...\n", +" S(1:N-1)'*flipdim(delta_d(k+J-N+1:k+J-1),1)+S(N)*d(k+J-N);\n", +" //Predicted value of the process\n", +" J=0;\n", +" yhat(k+J)=S(1:N-1)'*flipdim(delta_u(k+J-N+1:k+J-1),1)+S(N)*u(k+J-N);\n", +" //Corrected prediction for unforced case\n", +" ydottilde(:,k+1)=ydot(:,k+1)+ones(P1,1)*(y(k)-yhat(k));\n", +" \n", +" //Predicted unforced error Eqn20-52\n", +" edot(:,k+1)=yr(:,k+1)-ydottilde(:,k+1);\n", +" \n", +" //Control move\n", +" delta_u(k)=Kc1*edot(:,k+1);\n", +" u(k)=u(k-1)+delta_u(k);\n", +" \n", +"end\n", +"subplot(1,2,1);\n", +"plot(t,y,'black-');\n", +"subplot(1,2,2);\n", +"plot2d2(t,u);\n", +"//=============Part(b) MPC-2========//\n", +"n=T/Ts*2+1; //no. of discrete points in our domain of analysis\n", +"//Input initialization T is the Time for simulation\n", +"//We take a few u values in negative time to facilitate \n", +"//usage of step response model\n", +"u=zeros(n,1);\n", +"d=zeros(n,1);\n", +"delta_u=zeros(n,1);\n", +"//delta_u(101+2)=1; //Step change at t=2 min\n", +"u=cumsum(delta_u);\n", +"delta_d=zeros(n,1);\n", +"//delta_d(101+8)=0.15;//disturbance t=8 min\n", +"d=cumsum(delta_d);\n", +"y=zeros(1,n);//Actual values\n", +"yhat=zeros(1,n); //predicted value\n", +"ydot=zeros(P2,n); //Unforced predictions\n", +"ydottilde=zeros(P2,n); //Corrected unforced predictions\n", +"yr=ones(P2,n);//reference trajectory(same as setpoint)\n", +"edot=zeros(P2,n);//predicted unforced error\n", +"t=-(n-1)/2:Ts:(n-1)/2;\n", +"for k=(n-1)/2+1:n-P2\n", +" \n", +" //Unforced predictions\n", +" for J=1:P2\n", +" ydot(J,k+1)=S(J+1:N-1)'*flipdim(delta_u(k+J-N+1:k-1),1)+S(N)*u(k+J-N);\n", +" end\n", +" //Actual values of the process\n", +" J=0;\n", +" y(k+J)=S(1:N-1)'*flipdim(delta_u(k+J-N+1:k+J-1),1)+S(N)*u(k+J-N)+...\n", +" S(1:N-1)'*flipdim(delta_d(k+J-N+1:k+J-1),1)+S(N)*d(k+J-N);\n", +" //Predicted value of the process\n", +" J=0;\n", +" yhat(k+J)=S(1:N-1)'*flipdim(delta_u(k+J-N+1:k+J-1),1)+S(N)*u(k+J-N);\n", +" //Corrected prediction for unforced case\n", +" ydottilde(:,k+1)=ydot(:,k+1)+ones(P2,1)*(y(k)-yhat(k));\n", +" \n", +" //Predicted unforced error Eqn20-52\n", +" edot(:,k+1)=yr(:,k+1)-ydottilde(:,k+1);\n", +" \n", +" //Control move\n", +" delta_u(k)=Kc2(1,:)*edot(:,k+1);\n", +" u(k)=u(k-1)+delta_u(k);\n", +" \n", +"end\n", +"subplot(1,2,1);\n", +"plot(t,y,'-.');\n", +"set(gca(),'data_bounds',[0 60 0 1.25]); //putting bounds on display\n", +"l=legend('MPC(P=3,M=1)','MPC(P=4,M=2)',position=4);\n", +"xtitle('Process Output','Time(min)','$y$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;\n", +"subplot(1,2,2);\n", +"plot2d2(t,u,style=2);\n", +"//=============Part(b) PID========//\n", +"n=T/Ts*2+1; //no. of discrete points in our domain of analysis\n", +"//Input initialization T is the Time for simulation\n", +"//We take a few u values in negative time to facilitate \n", +"//usage of step response model\n", +"u=zeros(n,1);\n", +"d=zeros(n,1);\n", +"delta_u=zeros(n,1);\n", +"delta_d=zeros(n,1);\n", +"//delta_d(101+8)=0.15;//disturbance t=8 min\n", +"d=cumsum(delta_d);\n", +"y=zeros(n,1);//Actual values\n", +"ysp=1;//setpoint\n", +"e=zeros(n,1);//error\n", +"delta_e=zeros(n,1);//error\n", +"t=-(n-1)/2:Ts:(n-1)/2;\n", +"//PID settings\n", +"Kc=2.27;taui=16.6;tauD=1.49;\n", +"for k=(n-1)/2+1:n-1\n", +" //Actual values of the process\n", +" J=0;\n", +" y(k+J)=S(1:N-1)'*flipdim(delta_u(k+J-N+1:k+J-1),1)+S(N)*u(k+J-N)+...\n", +" S(1:N-1)'*flipdim(delta_d(k+J-N+1:k+J-1),1)+S(N)*d(k+J-N);\n", +" //error\n", +" e(k)=ysp-y(k);\n", +" delta_e(k)=e(k)-e(k-1);\n", +" \n", +" //Controller move----Digital PID----Eqn 7-28 Pg 136 (Velocity form)\n", +" u(k)=u(k-1)+Kc*([delta_e(k,1)+e(k,1)*Ts/taui+tauD/Ts*(e(k)-2*e(k-1)+e(k-2))]);\n", +" delta_u(k)=u(k)-u(k-1); \n", +"end\n", +"subplot(1,2,1);\n", +"plot(t,y,'red--');\n", +"set(gca(),'data_bounds',[0 60 0 1.25]); //putting bounds on display\n", +"l=legend('MPC (P=3,M=1)','MPC (P=4,M=2)','PID controller',position=4);\n", +"xtitle('Process Output','Time(min)','$y$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;\n", +"subplot(1,2,2);\n", +"plot2d2(t,u,style=5);\n", +"set(gca(),'data_bounds',[0 30 -100 100]); //putting bounds on display\n", +"l=legend('MPC (P=3,M=1)','MPC (P=4,M=2)','PID controller',position=4);\n", +"xtitle('Controller Output','Time(min)','$u$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;\n", +"//=============Part(c)========//\n", +"//=============Part(c)========//\n", +"//=============Part(c)========//\n", +"//=============Part(c)========//\n", +"//=============Part(c)========//\n", +"//=============Part(c)========//\n", +"//=============Part(c)========//\n", +"//=============Part(c) MPC-1========//\n", +"n=T/Ts*2+1; //no. of discrete points in our domain of analysis\n", +"//Input initialization T is the Time for simulation\n", +"//We take a few u values in negative time to facilitate \n", +"//usage of step response model\n", +"u=zeros(n,1);\n", +"d=zeros(n,1);\n", +"delta_u=zeros(n,1);\n", +"u=cumsum(delta_u);\n", +"delta_d=zeros(n,1);\n", +"delta_d((n-1)/2+1)=1;//disturbance t=0 min\n", +"d=cumsum(delta_d);\n", +"y=zeros(1,n);//Actual values\n", +"yhat=zeros(1,n); //predicted value\n", +"ydot=zeros(P1,n); //Unforced predictions\n", +"ydottilde=zeros(P1,n); //Corrected unforced predictions\n", +"yr=zeros(P1,n);//reference trajectory(same as setpoint)\n", +"edot=zeros(P1,n);//predicted unforced error\n", +"t=-(n-1)/2:Ts:(n-1)/2;\n", +"for k=(n-1)/2+1:n-P1\n", +" \n", +" //Unforced predictions\n", +" for J=1:P1\n", +" ydot(J,k+1)=S(J+1:N-1)'*flipdim(delta_u(k+J-N+1:k-1),1)+S(N)*u(k+J-N);\n", +" end\n", +" //Actual values of the process\n", +" J=0;\n", +" y(k+J)=S(1:N-1)'*flipdim(delta_u(k+J-N+1:k+J-1),1)+S(N)*u(k+J-N)+...\n", +" S(1:N-1)'*flipdim(delta_d(k+J-N+1:k+J-1),1)+S(N)*d(k+J-N);\n", +" //Predicted value of the process\n", +" J=0;\n", +" yhat(k+J)=S(1:N-1)'*flipdim(delta_u(k+J-N+1:k+J-1),1)+S(N)*u(k+J-N);\n", +" //Corrected prediction for unforced case\n", +" ydottilde(:,k+1)=ydot(:,k+1)+ones(P1,1)*(y(k)-yhat(k));\n", +" \n", +" //Predicted unforced error Eqn20-52\n", +" edot(:,k+1)=yr(:,k+1)-ydottilde(:,k+1);\n", +" \n", +" //Control move\n", +" delta_u(k)=Kc1*edot(:,k+1);\n", +" u(k)=u(k-1)+delta_u(k);\n", +" \n", +"end\n", +"scf();\n", +"subplot(1,2,1);\n", +"plot(t,y,'black-');\n", +"subplot(1,2,2);\n", +"plot2d2(t,u);\n", +"//=============Part(c) MPC-2========//\n", +"n=T/Ts*2+1; //no. of discrete points in our domain of analysis\n", +"//Input initialization T is the Time for simulation\n", +"//We take a few u values in negative time to facilitate \n", +"//usage of step response model\n", +"u=zeros(n,1);\n", +"d=zeros(n,1);\n", +"delta_u=zeros(n,1);\n", +"u=cumsum(delta_u);\n", +"delta_d=zeros(n,1);\n", +"delta_d((n-1)/2+1)=1;//disturbance t=0 min\n", +"d=cumsum(delta_d);\n", +"y=zeros(1,n);//Actual values\n", +"yhat=zeros(1,n); //predicted value\n", +"ydot=zeros(P2,n); //Unforced predictions\n", +"ydottilde=zeros(P2,n); //Corrected unforced predictions\n", +"yr=zeros(P2,n);//reference trajectory(same as setpoint)\n", +"edot=zeros(P2,n);//predicted unforced error\n", +"t=-(n-1)/2:Ts:(n-1)/2;\n", +"for k=(n-1)/2+1:n-P2\n", +" \n", +" //Unforced predictions\n", +" for J=1:P2\n", +" ydot(J,k+1)=S(J+1:N-1)'*flipdim(delta_u(k+J-N+1:k-1),1)+S(N)*u(k+J-N);\n", +" end\n", +" //Actual values of the process\n", +" J=0;\n", +" y(k+J)=S(1:N-1)'*flipdim(delta_u(k+J-N+1:k+J-1),1)+S(N)*u(k+J-N)+...\n", +" S(1:N-1)'*flipdim(delta_d(k+J-N+1:k+J-1),1)+S(N)*d(k+J-N);\n", +" //Predicted value of the process\n", +" J=0;\n", +" yhat(k+J)=S(1:N-1)'*flipdim(delta_u(k+J-N+1:k+J-1),1)+S(N)*u(k+J-N);\n", +" //Corrected prediction for unforced case\n", +" ydottilde(:,k+1)=ydot(:,k+1)+ones(P2,1)*(y(k)-yhat(k));\n", +" \n", +" //Predicted unforced error Eqn20-52\n", +" edot(:,k+1)=yr(:,k+1)-ydottilde(:,k+1);\n", +" \n", +" //Control move\n", +" delta_u(k)=Kc2(1,:)*edot(:,k+1);\n", +" u(k)=u(k-1)+delta_u(k);\n", +" \n", +"end\n", +"subplot(1,2,1);\n", +"plot(t,y,'-.');\n", +"subplot(1,2,2);\n", +"plot2d2(t,u,style=2);\n", +"//=============Part(C) PID========//\n", +"n=T/Ts*2+1; //no. of discrete points in our domain of analysis\n", +"//Input initialization T is the Time for simulation\n", +"//We take a few u values in negative time to facilitate \n", +"//usage of step response model\n", +"u=zeros(n,1);\n", +"d=zeros(n,1);\n", +"delta_u=zeros(n,1);\n", +"u=cumsum(delta_u);\n", +"delta_d=zeros(n,1);\n", +"delta_d((n-1)/2+1)=1;//disturbance t=0 min\n", +"d=cumsum(delta_d);\n", +"y=zeros(n,1);//Actual values\n", +"ysp=0;//setpoint\n", +"e=zeros(n,1);//error\n", +"delta_e=zeros(n,1);//error\n", +"t=-(n-1)/2:Ts:(n-1)/2;\n", +"//PID settings\n", +"Kc=3.52;taui=6.98;tauD=1.73;\n", +"for k=(n-1)/2+1:n-1\n", +" //Actual values of the process\n", +" J=0;\n", +" y(k+J)=S(1:N-1)'*flipdim(delta_u(k+J-N+1:k+J-1),1)+S(N)*u(k+J-N)+...\n", +" S(1:N-1)'*flipdim(delta_d(k+J-N+1:k+J-1),1)+S(N)*d(k+J-N);\n", +" //error\n", +" e(k)=ysp-y(k);\n", +" delta_e(k)=e(k)-e(k-1);\n", +" \n", +" //Controller move----Digital PID----Eqn 7-28 Pg 136 (Velocity form)\n", +" u(k)=u(k-1)+Kc*([delta_e(k,1)+e(k,1)*Ts/taui+tauD/Ts*(e(k)-2*e(k-1)+e(k-2))]);\n", +" delta_u(k)=u(k)-u(k-1); \n", +"end\n", +"subplot(1,2,1);\n", +"plot(t,y,'red--');\n", +"set(gca(),'data_bounds',[0 60 -0.1 0.3]); //putting bounds on display\n", +"l=legend('MPC (P=3,M=1)','MPC (P=4,M=2)','PID controller',position=1);\n", +"xtitle('Process Output','Time(min)','$y$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;\n", +"subplot(1,2,2);\n", +"plot2d2(t,u,style=5);\n", +"set(gca(),'data_bounds',[0 30 -1.5 0]); //putting bounds on display\n", +"l=legend('MPC (P=3,M=1)','MPC (P=4,M=2)','PID controller',position=1);\n", +"xtitle('Controller Output','Time(min)','$u$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/21-Process_Monitoring.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/21-Process_Monitoring.ipynb new file mode 100644 index 0000000..d3b7cb7 --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/21-Process_Monitoring.ipynb @@ -0,0 +1,389 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 21: Process Monitoring" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 21.2: Semiconductor_processing_control_charts.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 21.2\n", +"disp('Example 21.2')\n", +"//data\n", +"x=[209.6 207.6 211.1 \n", +" 183.5 193.1 202.4 \n", +" 190.1 206.8 201.6 \n", +" 206.9 189.3 204.1 \n", +" 260. 209. 212.2 \n", +" 193.9 178.8 214.5 \n", +" 206.9 202.8 189.7 \n", +" 200.2 192.7 202.1 \n", +" 210.6 192.3 205.9 \n", +" 186.6 201.5 197.4 \n", +" 204.8 196.6 225. \n", +" 183.7 209.7 208.6 \n", +" 185.6 198.9 191.5 \n", +" 202.9 210.1 208.1 \n", +" 198.6 195.2 150. \n", +" 188.7 200.7 207.6 \n", +" 197.1 204. 182.9 \n", +" 194.2 211.2 215.4 \n", +" 191. 206.2 183.9 \n", +" 202.5 197.1 211.1 \n", +" 185.1 186.3 188.9 \n", +" 203.1 193.1 203.9 \n", +" 179.7 203.3 209.7 \n", +" 205.3 190. 208.2 \n", +" 203.4 202.9 200.4 ]\n", +"//Original Limits\n", +"n=3;\n", +"xbar=sum(x,2)/n; //mean calculation\n", +"s=sqrt(1/(n-1)*sum((x-repmat(xbar,1,3)).^2,2)); //standard deviation calculation\n", +"p=length(xbar);//no. of subgroups\n", +"xbarbar=sum(xbar,1)/p;\n", +"sbar=sum(s,1)/p;\n", +"c4=0.8862;B3=0;B4=2.568;c=3;\n", +"sigma=1/c4*sbar/sqrt(n);\n", +"//original limits\n", +"UCL_x=xbarbar+c*sigma;//Eqn21-9\n", +"LCL_x=xbarbar-c*sigma;//Eqn 21-10\n", +"UCL_s=B4*sbar;//Eqn21-14\n", +"LCL_s=B3*sbar;//Eqn21-15\n", +"//Modified Limits\n", +"x_mod=x;\n", +"x_mod([5 15],:)=[];\n", +"n=3;\n", +"xbar_mod=sum(x_mod,2)/n; //mean calculation\n", +"s_mod=sqrt(1/(n-1)*sum((x_mod-repmat(xbar_mod,1,3)).^2,2)); //standard deviation calculation\n", +"p_mod=length(xbar_mod);//no. of subgroups\n", +"xbarbar_mod=sum(xbar_mod,1)/p_mod;\n", +"sbar_mod=sum(s_mod,1)/p_mod;\n", +"c4=0.8862;B3=0;B4=2.568;c=3;\n", +"sigma_mod=1/c4*sbar_mod/sqrt(n);\n", +"//modified limits\n", +"UCL_x_mod=xbarbar_mod+c*sigma_mod;//Eqn21-9\n", +"LCL_x_mod=xbarbar_mod-c*sigma_mod;//Eqn 21-10\n", +"UCL_s_mod=B4*sbar_mod;//Eqn21-14\n", +"LCL_s_mod=B3*sbar_mod;//Eqn21-15\n", +"mprintf('\n Original Limits Modified Limits')\n", +"mprintf('\n xbar Chart Control Limits ')\n", +"mprintf('\n UCL %f %f',UCL_x,UCL_x_mod)\n", +"mprintf('\n LCL %f %f',LCL_x,LCL_x_mod)\n", +"mprintf('\n s Chart Control Limits ')\n", +"mprintf('\n UCL %f %f',UCL_s,UCL_s_mod)\n", +"mprintf('\n LCL %f %f',LCL_s,LCL_s_mod)\n", +"subplot(2,1,1);\n", +"plot2d(repmat(UCL_x,1,p));\n", +"plot(repmat(UCL_x_mod,1,p),'--');\n", +"plot2d(repmat(LCL_x,1,p));\n", +"plot(repmat(LCL_x_mod,1,p),'--');\n", +"plot2d(xbar,style=-1,rect=[0,160,30,260])\n", +"xtitle('Example 21.2','Sampling number','$\text{Thickness}\ (\AA)$')\n", +"legend('Original Limits','Modified Limits')\n", +"subplot(2,1,2);\n", +"plot2d(repmat(UCL_s,1,p));\n", +"plot2d(repmat(LCL_s,1,p));\n", +"plot(repmat(UCL_s_mod,1,p),'--');\n", +"plot(repmat(LCL_s_mod,1,p),'--');\n", +"plot2d(s,style=-1,rect=[0,0,30,30])\n", +"xtitle('','Sampling number','$\text{Standard Deviation}\ (\AA)$')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 21.3: Shewart_CUSUM_EWMA_compariso.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 21.3\n", +"disp('Example 21.3')\n", +"T=70;sigma=3;p=100;//p is the no. of samples\n", +"x=grand(p,1, 'nor', T, sigma);\n", +"delta=0.5*sigma;\n", +"x(11:$)=x(11:$)+delta;\n", +"//Limits for Shewart charts\n", +"UCL_1=T+sigma*3;\n", +"LCL_1=T-sigma*3;\n", +"subplot(3,1,1);\n", +"plot2d(repmat(UCL_1,1,p));\n", +"plot2d(repmat(LCL_1,1,p));\n", +"plot2d(x,style=1,rect=[0,60,100,80])\n", +"xtitle('Example 21.3','Sampling number','Strength (MPa)')\n", +"//CUSUM\n", +"Cplus=zeros(100,1);Cminus=zeros(100,1);\n", +"K=0.5*sigma;H=5*sigma;\n", +"UCL_2=H;\n", +"for k=2:100\n", +" Cplus(k)=max(0,x(k)-(T+K)+Cplus(k-1));\n", +" Cminus(k)=max(0,(T-K)-x(k)+Cminus(k-1));\n", +" if Cplus(k-1)>H then \n", +" Cplus(k)=0;\n", +" end\n", +" if Cminus(k-1)>H then \n", +" Cminus(k)=0;\n", +" end\n", +" \n", +"end\n", +"subplot(3,1,2);\n", +"plot2d(Cplus,style=1,rect=[0,0,100,20]);\n", +"plot2d(Cminus,style=2,rect=[0,0,100,20]);\n", +"plot2d(repmat(UCL_2,1,p));\n", +"xtitle('','Sampling number','CUSUM')\n", +"legend('Cplus','Cminus')\n", +"//EWMA\n", +"lamda=0.25;\n", +"z=x;\n", +"for k=2:100\n", +" z(k)=lamda*x(k)+(1-lamda)*z(k-1);\n", +"end\n", +"UCL_3=T+3*sigma*sqrt(lamda/(2-lamda));\n", +"LCL_3=T-3*sigma*sqrt(lamda/(2-lamda));\n", +"subplot(3,1,3);\n", +"plot2d(repmat(UCL_3,1,p));\n", +"plot2d(repmat(LCL_3,1,p));\n", +"plot2d(z,style=1,rect=[0,65,100,75])\n", +"xtitle('','Sampling number','EWMA')\n", +"mprintf('The charts in the example and in the book differ due\n...\n", +"a different realization of data everytime the code is run\n...\n", +"due to the grand command. If we had the exact data as that given\n...\n", +"in the book our charts would have matched.')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 21.4: Process_Capability_Indices.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 21.4\n", +"disp('Example 21.4')\n", +"xbar=199.5;//Note that this is the correct value and not 199\n", +"sbar=8.83;\n", +"USL=235;//Note that this is diff from UCL\n", +"LSL=185;\n", +"c4=0.8862;\n", +"n=3;\n", +"sigma=5.75;\n", +"sigma_x=sbar/c4/sqrt(n);\n", +"mprintf('\nValue of sigma_x=%f',sigma_x);\n", +"Cp=(USL-LSL)/6/sigma;\n", +"Cpk=min(xbar-LSL,USL-xbar)/3/sigma;\n", +"mprintf('\nCp=%f and Cpk=%f',Cp,Cpk)" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 21.5: Effluent_Stream_from_wastewater_treatment.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 21.5\n", +"disp('Example 21.5')\n", +"//data\n", +"x=[ 17.7 1380. \n", +" 23.6 1458. \n", +" 13.2 1322. \n", +" 25.2 1448. \n", +" 13.1 1334. \n", +" 27.8 1485. \n", +" 29.8 1503. \n", +" 9. 1540. \n", +" 14.3 1341. \n", +" 26. 1448. \n", +" 23.2 1426. \n", +" 22.8 1417. \n", +" 20.4 1384. \n", +" 17.5 1380. \n", +" 18.4 1396. \n", +" 16.8 1345. \n", +" 13.8 1349. \n", +" 19.4 1398. \n", +" 24.7 1426. \n", +" 16.8 1361. \n", +" 14.9 1347. \n", +" 27.6 1476. \n", +" 26.1 1454. \n", +" 20. 1393. \n", +" 22.9 1427. \n", +" 22.4 1431. \n", +" 19.6 1405. \n", +" 31.5 1521. \n", +" 19.9 1409. \n", +" 20.3 1392.]; \n", +"T=mean(x,'r');\n", +"s=sqrt(variance(x,'r'));\n", +"UCL=T+3*s;\n", +"LCL=T-3*s;\n", +"p=size(x,1)\n", +"subplot(2,1,1);\n", +"plot2d(repmat(UCL(1),1,p));\n", +"plot2d(repmat(LCL(1),1,p));\n", +"plot2d(repmat(T(1),1,p));\n", +"plot2d(x(:,1),style=-1,rect=[0,0,32,40])\n", +"xtitle('Example 21.4','Sampling number','BOD (mg/L)')\n", +"subplot(2,1,2);\n", +"plot2d(repmat(UCL(2),1,p));\n", +"plot2d(repmat(LCL(2),1,p));\n", +"plot2d(repmat(T(2),1,p));\n", +"plot2d(x(:,2),style=-1,rect=[0,1200,32,1600])\n", +"xtitle('','Sampling number','Solids (mg/L)')\n", +"//subplot(3,1,3);\n", +"scf()\n", +"plot2d(x(8,1),x(8,2),style=-3,rect=[0,1200,40,1600])\n", +"plot2d(x(:,1),x(:,2),style=-1,rect=[0,1200,40,1600])\n", +"legend('Sample #8',position=4)\n", +"xtitle('','BOD (mg/L)','Solids (mg/L)')\n", +"mprintf('\nThe confidence interval for third case is not drawn\n...\n", +" because it is beyond the scope of this book')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 21.6: Hotellings_T_square_statistic.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 21.6\n", +"disp('Example 21.6')\n", +"//data\n", +"x=[ 17.7 1380. \n", +" 23.6 1458. \n", +" 13.2 1322. \n", +" 25.2 1448. \n", +" 13.1 1334. \n", +" 27.8 1485. \n", +" 29.8 1503. \n", +" 9. 1540. \n", +" 14.3 1341. \n", +" 26. 1448. \n", +" 23.2 1426. \n", +" 22.8 1417. \n", +" 20.4 1384. \n", +" 17.5 1380. \n", +" 18.4 1396. \n", +" 16.8 1345. \n", +" 13.8 1349. \n", +" 19.4 1398. \n", +" 24.7 1426. \n", +" 16.8 1361. \n", +" 14.9 1347. \n", +" 27.6 1476. \n", +" 26.1 1454. \n", +" 20. 1393. \n", +" 22.9 1427. \n", +" 22.4 1431. \n", +" 19.6 1405. \n", +" 31.5 1521. \n", +" 19.9 1409. \n", +" 20.3 1392.]; \n", +" \n", +" \n", +"n=1;\n", +"N=size(x,1);\n", +"T=mean(x,'r');\n", +"//For our example n=1 because each measurement is a subgroup\n", +"S=mvvacov(x);\n", +"//Note that mvvacov calculates covariance with denominator N, while\n", +"//variance caluclates with denominator N-1, hence diagonal entry of mvvacov does not\n", +"//match with variance calculated manually for each vector\n", +"//As per wikipedia the book is wrong and for covariance matrix we should\n", +"//use N-1 but here we follow the book\n", +"Tsquare=zeros(N,1);\n", +"for k=1:N\n", +" Tsquare(k)=n*(x(k,:)-T)*inv(S)*(x(k,:)-T)';\n", +"end\n", +"UCL=11.63;\n", +"plot(repmat(UCL,1,N),color='black');\n", +"plot(Tsquare,'+')\n", +"legend('UCL 99% confidence limit')\n", +"xtitle('Example 21.6','Sample number','$T^2$')" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/22-Biosystems_Control_Design.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/22-Biosystems_Control_Design.ipynb new file mode 100644 index 0000000..2270484 --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/22-Biosystems_Control_Design.ipynb @@ -0,0 +1,447 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 22: Biosystems Control Design" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 22.1: Fermentor.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 22.1\n", +"disp('Example 22.1')\n", +"//Parameters\n", +"Yxs=0.4;B=0.2;Pm=50;Ki=22;\n", +"a=2.2;mu_m=0.48;Km=1.2;Sf=20;\n", +"//ODE model\n", +"function ydot=model(t,y,D)\n", +" X=y(1);S=y(2);P=y(3);\n", +" \n", +" Xdot=-D*X+mu(S,P)*X;\n", +" Sdot=D*(Sf-S)-1/Yxs*mu(S,P)*X;\n", +" Pdot=-D*P+[a*mu(S,P)+B]*X\n", +" \n", +" ydot=[Xdot Sdot Pdot]';\n", +"endfunction\n", +"//Rate law\n", +"function mu=mu(S,P)\n", +" mu=mu_m*(1-P/Pm)*S/(Km+S+S^2/Ki); \n", +"endfunction\n", +"t=0:0.1:100;t0=0;\n", +"y0=[6 5 19.14]';//Initial stable condition\n", +"D=0.202*1.1;//10% increase\n", +"y_up = ode(y0, t0, t, list(model,D))\n", +"D=0.202*0.9;//10% decrease\n", +"y_down = ode(y0, t0, t, list(model,D))\n", +"subplot(2,1,1);\n", +"plot(t,y_up(1,:),color='red');\n", +"plot(t,y_down(1,:));\n", +"xtitle('$D=0.202\ h^{-1}$','Time(h)','Biomass (g/L)')\n", +"legend('Dilution +10%','Dilution -10%')\n", +"subplot(2,1,2);\n", +"t=0:0.1:100;t0=0;\n", +"y0=[6 5 44.05]';//Initial stable condition\n", +"D=0.0389*1.1;//10% increase\n", +"y_up = ode(y0, t0, t, list(model,D))\n", +"D=0.0389*0.9;//10% decrease\n", +"y_down = ode(y0, t0, t, list(model,D))\n", +"plot(t,y_up(1,:),color='red');\n", +"plot(t,y_down(1,:))\n", +"xtitle('$D=0.0389\ h^{-1}$','Time(h)','Biomass (g/L)');\n", +"legend('Dilution +10%','Dilution -10%')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 22.2: Granulation_process_control.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 22.2\n", +"disp('Example 22.2')\n", +"//Author: Dhruv Gupta....Aug 4, 2013\n", +"//<dgupta6@wisc.edu>\n", +"K=[0.2 0.58 0.35;0.25 1.10 1.3;0.3 0.7 1.2];\n", +"tau=[2 2 2;3 3 3;4 4 4];\n", +"s=%s;\n", +"G=K./(1+tau*s);\n", +"RGA=K.*inv(K');\n", +"disp(RGA,'RGA=')\n", +"//IMC based tuning\n", +"tauC=5;\n", +"Kc=diag(tau/tauC./K);\n", +"mprintf('\n\nThe tauI given in book are wrong\n...\n", +"refer to Table 11.1 for correct formula\n\n')\n", +"tauI=diag(tau)+1; \n", +"mprintf('\nWe still however use the ones given in book\n');\n", +"disp(Kc,'Kc=')\n", +"disp(tauI,'tauI=')\n", +"//Refer to Eqns 15-23 and 15-24\n", +"Gc=Kc.*(1+(1)./tauI/s);\n", +"//For the sake of brevity we write Gstar as G\n", +"//We will account for delays in the for loop that we will write\n", +"//Refer to Figure 15.9 Page 295 for details of Smith Predictor\n", +"//====Making step response models of the continuos transfer functions====//\n", +"Ts=0.1;//Sampling time ie delta_T\n", +"delay=3/Ts;\n", +"N=150/Ts;//Model Order\n", +"s=%s;\n", +"G=syslin('c',diag(matrix(G,1,9)));//Transfer function\n", +"t=0:Ts:N*Ts;\n", +"u_sim=ones(9,length(t));\n", +"//u_sim(:,1:4)=zeros(9,4); //input delay to account for 3 min delay in G\n", +"S=csim(u_sim,t,G)';//generating step response model for real plant\n", +"//plot(t,S);\n", +"S(1,:)=[];\n", +"//Now we have these step response models for each of the transfer functions\n", +"//[S1 S4 S7\n", +"//S2 S5 S8\n", +"//S3 S6 S9]\n", +"T=150+delay;//Simulation Run Time in minutes(we add delay because our for loop runs till n-delay)\n", +"n=T/Ts*2+1; //no. of discrete points in our domain of analysis\n", +"//Input initialization T is the Time for simulation\n", +"//========Set point as 10=============//\n", +"//p is the controller output\n", +"p=zeros(n,3);\n", +"delta_p=zeros(n,3);\n", +"ytilde=zeros(n,3); //Prediction of Smith Fig 15.9\n", +"e=zeros(n,3); //corrections\n", +"edash=zeros(n,3);\n", +"delta_edash=zeros(n,3);\n", +"ysp=zeros(n,3);\n", +"ysp((n-1)/2+1:n,1)=10*ones(n-((n-1)/2+1)+1,1);\n", +"t=-(n-1)/2*Ts:Ts:(n-1)/2*Ts;\n", +"y=zeros(n,3);\n", +"for k=(n-1)/2+1:n-delay\n", +" \n", +" //Error e\n", +" e(k,:)=ysp(k-1,:)-y(k-1,:);\n", +" \n", +" //Error edash\n", +" edash(k,:)=e(k-1,:)-ytilde(k-1,:)+ytilde(k-1-delay,:);\n", +" //Edash=E-(Y1-Y2)...where Y2 is delayed Y1\n", +" delta_edash(k,:)=edash(k,:)-edash(k-1,:);\n", +" //Controller calculation----Digital PID----Eqn 7-28 Pg 136 (Velocity form)\n", +" p(k,:)=p(k-1,:)+[delta_edash(k,:)+edash(k,:)*diag(Ts./tauI)]*diag(Kc);\n", +" \n", +" //Limits on manipulated variables\n", +" p(k,:)=min([(345-180)*ones(1,3);p(k,:)],'r');\n", +" p(k,:)=max([(105-180)*ones(1,3);p(k,:)],'r');\n", +" \n", +" delta_p(k,:)=p(k,:)-p(k-1,:);\n", +" \n", +" \n", +" //Prediction\n", +" ytilde(k,1)=[S(1:N-1,1);S(1:N-1,4);S(1:N-1,7)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);flipdim(delta_p(k-N+1:k-1,3),1)]...\n", +" +[S(N,1) S(N,4) S(N,7)]*[p(k-N,1);p(k-N,2);p(k-N,3)];\n", +" ytilde(k,2)=[S(1:N-1,2);S(1:N-1,5);S(1:N-1,8)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);flipdim(delta_p(k-N+1:k-1,3),1)]...\n", +" +[S(N,2) S(N,5) S(N,8)]*[p(k-N,1);p(k-N,2);p(k-N,3)];\n", +" ytilde(k,3)=[S(1:N-1,3);S(1:N-1,6);S(1:N-1,9)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);flipdim(delta_p(k-N+1:k-1,3),1)]...\n", +" +[S(N,3) S(N,6) S(N,9)]*[p(k-N,1);p(k-N,2);p(k-N,3)];\n", +" \n", +" //Output\n", +" y(k+delay,1)=[S(1:N-1,1);S(1:N-1,4);S(1:N-1,7)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);flipdim(delta_p(k-N+1:k-1,3),1)]...\n", +" +[S(N,1) S(N,4) S(N,7)]*[p(k-N,1);p(k-N,2);p(k-N,3)];\n", +" y(k+delay,2)=[S(1:N-1,2);S(1:N-1,5);S(1:N-1,8)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);flipdim(delta_p(k-N+1:k-1,3),1)]...\n", +" +[S(N,2) S(N,5) S(N,8)]*[p(k-N,1);p(k-N,2);p(k-N,3)];\n", +" y(k+delay,3)=[S(1:N-1,3);S(1:N-1,6);S(1:N-1,9)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);flipdim(delta_p(k-N+1:k-1,3),1)]...\n", +" +[S(N,3) S(N,6) S(N,9)]*[p(k-N,1);p(k-N,2);p(k-N,3)];\n", +"end\n", +"subplot(2,2,1);\n", +"plot(t',y(:,1),'--',t',y(:,2),':',t',y(:,3),'-.',t',ysp(:,1),'-');\n", +"set(gca(),'data_bounds',[0 150 -4 12]); //putting bounds on display\n", +"l=legend('$y1$','$y2$','$y3$',position=1);\n", +"l.font_size=5;\n", +"xtitle('','Time(min)','$y$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;\n", +"subplot(2,2,2);\n", +"plot(t',p(:,1),'--',t',p(:,2),':',t',p(:,3),'-.');\n", +"set(gca(),'data_bounds',[-1 150 -40 100]); //putting bounds on display\n", +"l=legend('$p1$','$p2$','$p3$',position=1);\n", +"l.font_size=5;\n", +"xtitle('','Time(min)','$p$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;\n", +"mprintf('Note that there is no overshoot around time=25 mins \n...\n", +"which is in contrast to what is shown in book')\n", +"//========Now for set point as 50=============//\n", +"//p is the controller output\n", +"p=zeros(n,3);\n", +"delta_p=zeros(n,3);\n", +"ytilde=zeros(n,3); //Prediction of Smith Fig 15.9\n", +"e=zeros(n,3); //corrections\n", +"edash=zeros(n,3);\n", +"delta_edash=zeros(n,3);\n", +"ysp=zeros(n,3);\n", +"ysp((n-1)/2+1:n,1)=50*ones(n-((n-1)/2+1)+1,1);\n", +"t=-(n-1)/2*Ts:Ts:(n-1)/2*Ts;\n", +"y=zeros(n,3);\n", +"for k=(n-1)/2+1:n-delay\n", +" \n", +" //Error e\n", +" e(k,:)=ysp(k-1,:)-y(k-1,:);\n", +" \n", +" //Error edash\n", +" edash(k,:)=e(k-1,:)-ytilde(k-1,:)+ytilde(k-1-delay,:);\n", +" //Edash=E-(Y1-Y2)...where Y2 is delayed Y1\n", +" delta_edash(k,:)=edash(k,:)-edash(k-1,:);\n", +" //Controller calculation----Digital PID----Eqn 7-28 Pg 136 (Velocity form)\n", +" p(k,:)=p(k-1,:)+[delta_edash(k,:)+edash(k,:)*diag(Ts./tauI)]*diag(Kc);\n", +" \n", +" //Limits on manipulated variables\n", +" p(k,:)=min([(345-180)*ones(1,3);p(k,:)],'r');\n", +" p(k,:)=max([(105-180)*ones(1,3);p(k,:)],'r');\n", +" \n", +" delta_p(k,:)=p(k,:)-p(k-1,:);\n", +" \n", +" \n", +" //Prediction\n", +" ytilde(k,1)=[S(1:N-1,1);S(1:N-1,4);S(1:N-1,7)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);flipdim(delta_p(k-N+1:k-1,3),1)]...\n", +" +[S(N,1) S(N,4) S(N,7)]*[p(k-N,1);p(k-N,2);p(k-N,3)];\n", +" ytilde(k,2)=[S(1:N-1,2);S(1:N-1,5);S(1:N-1,8)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);flipdim(delta_p(k-N+1:k-1,3),1)]...\n", +" +[S(N,2) S(N,5) S(N,8)]*[p(k-N,1);p(k-N,2);p(k-N,3)];\n", +" ytilde(k,3)=[S(1:N-1,3);S(1:N-1,6);S(1:N-1,9)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);flipdim(delta_p(k-N+1:k-1,3),1)]...\n", +" +[S(N,3) S(N,6) S(N,9)]*[p(k-N,1);p(k-N,2);p(k-N,3)];\n", +" \n", +" //Output\n", +" y(k+delay,1)=[S(1:N-1,1);S(1:N-1,4);S(1:N-1,7)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);flipdim(delta_p(k-N+1:k-1,3),1)]...\n", +" +[S(N,1) S(N,4) S(N,7)]*[p(k-N,1);p(k-N,2);p(k-N,3)];\n", +" y(k+delay,2)=[S(1:N-1,2);S(1:N-1,5);S(1:N-1,8)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);flipdim(delta_p(k-N+1:k-1,3),1)]...\n", +" +[S(N,2) S(N,5) S(N,8)]*[p(k-N,1);p(k-N,2);p(k-N,3)];\n", +" y(k+delay,3)=[S(1:N-1,3);S(1:N-1,6);S(1:N-1,9)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);flipdim(delta_p(k-N+1:k-1,3),1)]...\n", +" +[S(N,3) S(N,6) S(N,9)]*[p(k-N,1);p(k-N,2);p(k-N,3)];\n", +"end\n", +"subplot(2,2,3);\n", +"plot(t',y(:,1),'--',t',y(:,2),':',t',y(:,3),'-.',t',ysp(:,1),'-');\n", +"set(gca(),'data_bounds',[0 150 -10 60]); //putting bounds on display\n", +"l=legend('$y1$','$y2$','$y3$',position=1);\n", +"l.font_size=5;\n", +"xtitle('','Time(min)','$y$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;\n", +"subplot(2,2,4);\n", +"plot(t',p(:,1),'--',t',p(:,2),':',t',p(:,3),'-.');\n", +"set(gca(),'data_bounds',[-1 150 -100 200]); //putting bounds on display\n", +"l=legend('$p1$','$p2$','$p3$',position=1);\n", +"l.font_size=5;\n", +"xtitle('','Time(min)','$p$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 22.3: Type_1_Diabetes.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 22.3\n", +"disp('Example 22.3')\n", +"//Parameters\n", +"p1=0.028735;p2=0.028344;p3=5.035E-5;V1=12;n=0.0926;\n", +"Ib=15;//basal\n", +"Gb=81;\n", +"//Diet function\n", +"function D=D(t)\n", +" D=9*exp(-0.05*t);\n", +"endfunction\n", +"//ODE model\n", +"function ydot=model(t,y,U)\n", +" G=y(1);X=y(2);I=y(3);\n", +" Gdot=-p1*G-X*(G+Gb)+D(t);\n", +" Xdot=-p2*X+p3*I;\n", +" Idot=-n*(I+Ib)+U/V1;\n", +" ydot=[Gdot Xdot Idot]';\n", +"endfunction\n", +"t=0:0.1:400;t0=0;\n", +"y0=[0 0 0]';//G,X,I are deviation variables\n", +"U=0;\n", +"y = Gb+ode(y0, t0, t, list(model,U))\n", +"subplot(3,1,1);\n", +"plot(t,y(1,:));\n", +"xtitle('','Time(min)','Glucose (mg/L)')\n", +"legend('$U=0\ mU/min$')\n", +"U=15;\n", +"y =Gb+ ode(y0, t0, t, list(model,U))\n", +"subplot(3,1,2);\n", +"plot(t,y(1,:));\n", +"xtitle('','Time(min)','Glucose (mg/L)')\n", +"legend('$U=15\ mU/min$')\n", +"U=25;\n", +"y = Gb+ode(y0, t0, t, list(model,U))\n", +"subplot(3,1,3);\n", +"plot(t,y(1,:));\n", +"xtitle('','Time(min)','Glucose (mg/L)')\n", +"legend('$U=25\ mU/min$')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 22.4: Influence_of_drugs.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 22.4\n", +"disp('Example 22.4')\n", +"K=[-6 3;12 5];\n", +"tau=[0.67 2;0.67 5];\n", +"s=%s;\n", +"G=K./(1+tau*s);\n", +"delay75=0.75;\n", +"delay1=1;\n", +"RGA=K.*inv(K');\n", +"disp(RGA,'RGA=')\n", +"//IMC based tuning\n", +"tauC=[tau(1,1) tau(2,2)];\n", +"Kc=diag(tau./(repmat(tauC,2,1)+[delay75 delay1;delay75 delay1])./K);\n", +"tauI=diag(tau);\n", +"disp(Kc,'Kc=')\n", +"disp(tauI,'tauI=')\n", +"//====Making step response models of the continuos transfer functions====//\n", +"Ts=0.05;//Sampling time ie delta_T\n", +"delay75=0.75/Ts;\n", +"delay1=1/Ts;\n", +"N=30/Ts;//Model Order\n", +"s=%s;\n", +"G=syslin('c',diag(matrix(G,1,4)));//Transfer function\n", +"t=0:Ts:N*Ts;\n", +"u_sim=ones(4,length(t));\n", +"u_sim([1 2],1:(delay75))=zeros(2,delay75); //input delay to account for delay in SNP\n", +"u_sim([3 4],1:(delay1))=zeros(2,delay1); //input delay to account for delay in DPM\n", +"S=csim(u_sim,t,G)';//generating step response model for real plant\n", +"//plot(t,S);\n", +"S(1,:)=[];\n", +"//Now we have these step response models for each of the transfer functions\n", +"//[S1 S3\n", +"//S2 S4 \n", +"T=50;//Simulation Run Time in minutes\n", +"n=T/Ts*2+1; //no. of discrete points in our domain of analysis\n", +"//========Set point as -10=============//\n", +"//p is the controller output\n", +"p=zeros(n,2);\n", +"delta_p=zeros(n,2);\n", +"e=zeros(n,2); //errors=(ysp-y) on which PI acts\n", +"ysp=zeros(n,2);\n", +"ysp((n-1)/2+1:n,1)=-10*ones(n-((n-1)/2+1)+1,1);\n", +"t=-(n-1)/2*Ts:Ts:(n-1)/2*Ts;\n", +"y=zeros(n,2);\n", +"for k=(n-1)/2+1:n\n", +" \n", +" //Error e\n", +" e(k,:)=ysp(k-1,:)-y(k-1,:);\n", +" delta_e(k,:)=e(k,:)-e(k-1,:);\n", +" //Controller calculation----Digital PID----Eqn 7-28 Pg 136 (Velocity form)\n", +"// p(k,:)=p(k-1,:)+flipdim([delta_e(k,:)+e(k,:)*diag(Ts./tauI)]*diag(Kc),2);\n", +" p(k,:)=p(k-1,:)+([delta_e(k,:)+e(k,:)*diag(Ts./tauI)]*diag(Kc));\n", +" //1-1/2-2 pairing\n", +" \n", +" delta_p(k,:)=p(k,:)-p(k-1,:);\n", +" \n", +" //Output\n", +" y(k,1)=[S(1:N-1,1);S(1:N-1,3)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1)]...\n", +" +[S(N,1) S(N,3)]*[p(k-N,1);p(k-N,2)];\n", +" y(k,2)=[S(1:N-1,2);S(1:N-1,4)]'...\n", +" *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);]...\n", +" +[S(N,2) S(N,4)]*[p(k-N,1);p(k-N,2)];\n", +"end\n", +"plot(t',y(:,1),'--',t',y(:,2),':',t',ysp(:,1),'-');\n", +"set(gca(),'data_bounds',[0 50 -20 25]); //putting bounds on display\n", +"l=legend('Mean arterial pressure','Cardiac output',position=1);\n", +"//l.font_size=5;\n", +"xtitle('','Time(min)','$y$');\n", +"a=get('current_axes');\n", +"c=a.y_label;c.font_size=5;\n", +"mprintf('\nThere is more interaction in the variables\n...\n", +"than the book claims, hence a mismatch between the result\n...\n", +"and the book\n')" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/6-Development_of_Empirical_Models_from_Process_Data.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/6-Development_of_Empirical_Models_from_Process_Data.ipynb new file mode 100644 index 0000000..0111bc4 --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/6-Development_of_Empirical_Models_from_Process_Data.ipynb @@ -0,0 +1,341 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 6: Development of Empirical Models from Process Data" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 6.1: Gas_turbine_generator.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 6.1\n", +"disp('Example 6.1')\n", +"//Fuel flow rate appended with ones for intercept in regression\n", +"fuel=[1 2.3 2.9 4 4.9 5.8 6.5 7.7 8.4 9]; \n", +"X=[ones(1,10);fuel]'; \n", +"Y=[2 4.4 5.4 7.5 9.1 10.8 12.3 14.3 15.8 16.8]'; //Power generated\n", +"Bhat=inv(X'*X)*X'*Y;\n", +"mprintf('\n Linear model \n B1_hat=%f \n B2_hat=%f',Bhat')\n", +"//For better accuracy we can fit higher order model\n", +"X_new=[ones(1,10);fuel;fuel.^2]';\n", +"Bhat_new=inv(X_new'*X_new)*X_new'*Y;\n", +"mprintf('\n \n Quadratic model \n B1_hat=%f \n B2_hat=%f \n B3_hat=%f',Bhat_new')\n", +"Output_table=[fuel' Y X*Bhat X_new*Bhat_new];\n", +"//mprintf('\n Fuel Power Generated Linear Model Quadratic Model %f %f',Output_table(:,1),Output_table(:,2))\n", +"//disp(Output_table)\n", +"//Table 6.1\n", +"mprintf('\n \n Table 6.1 %s','')\n", +"mprintf('\n ui yi Linear Model Quadratic Model %s','')\n", +"mprintf('\n %f %f %f %15f',Output_table)\n", +"//Error calculations ----(This is not given in book-requires understanding of statistics)\n", +"Yhat=X*Bhat; //Predicted Y from regression variables\n", +"S_lin=(Y-Yhat)'*(Y-Yhat); //Sum of errors in Y for linear model --eqn 6.9\n", +"S_quad=(Y-X_new*Bhat_new)'*(Y-X_new*Bhat_new); //Errors in Y for quadratic model\n", +"mprintf('\n %25s%f %10s%f','S=',S_lin,'S=',S_quad)\n", +"n=length(fuel);\n", +"sigma=S_lin/(n-1)*(inv(X'*X));\n", +"bounds=(sigma.^0.5)/sqrt(n)*2.262;\n", +"mprintf('\n The errors in Bhats are not calculated because the procedure is not...\n", +"\n given in the solution of the example')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 6.2: Continuous_stirred_tank_reactor.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 6.2\n", +"disp('Example 6.2')\n", +"deltaw=5;//kg/min\n", +"deltaT=20;//deg C\n", +"K=deltaT/deltaw\n", +"tau=5//min\n", +"T=140+0.632*20;//152.6 deg C\n", +"s=%s;\n", +"G=4/(5*s+1); //G=T'(s)/W'(s)\n", +"mprintf('\n T(s)/W(s)=%s','')\n", +"disp(G)\n", +"t=0:0.01:25;\n", +"n=length(t);\n", +"w=5*ones(1,n);\n", +"yt=csim(w,t,G)+140;\n", +"wt=w+120;\n", +"subplot(2,1,2);\n", +"plot2d(t,yt,rect=[0,130,25,160]);\n", +"xtitle('','Time(min)','$T(^0C)$')\n", +"xgrid();\n", +"subplot(2,1,1);\n", +"plot2d(t,wt,rect=[0,120,25,126],style=2)\n", +"xtitle('Fig 6.4','Time(min)','$w\ (kg/min)$')\n", +"xgrid();" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 6.3: Off_gas_C02_concentration_response.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 6.3\n", +"disp('Example 6.3')\n", +"//Smith's method\n", +"t20=1.85;//min\n", +"t60=5;//min\n", +"ratio=t20/t60;\n", +"zeta=1.3;//Zeta obtained from Fig 6.7 page 109\n", +"tau=t60/2.8//Value 2.8 obtained from Fig 6.7\n", +"tau1=tau*zeta+tau*sqrt(zeta^2-1);\n", +"tau2=tau*zeta-tau*sqrt(zeta^2-1);\n", +"mprintf('From Smiths method \n tau1=%f min\n tau2=%f min \n',[tau1 tau2])\n", +"//Nonlinear regression\n", +"//This method cannot be directly used here because we do not have the data with us\n", +"//Had the data been given in tabular form we could have performed a regression\n", +"//Converting graphical data(Fig 7.8) printed in textbook to tabular form is not practical\n", +"//Towards the end of the program however a roundabout way has been implemented so\n", +"//that the user can learn the technique of non-linear optimization\n", +"s=%s;\n", +"G2=1/((tau1*s+1)*(tau2*s+1)) //Smith's method\n", +"G3=1/(4.60*s+1)//First order with time delay: Note that we cannot have exp(-0.7s) without symbolic toolbox so we use a roundaround trick for this\n", +"//Also note that we could have used Pade's approximation but that works well only for very small time delays\n", +"G1=1/((3.34*s+1)*(1.86*s+1)); //Nonlinear regression\n", +"Glist=syslin('c',[G1 G2 G3]') //Simply collating them together for further simulation\n", +"G=syslin('c',Glist);\n", +"t=0:0.1:20;\n", +"y=csim('step',t,G); \n", +"y(3,:)=[zeros(1,8) y(3,1:$-8)] //This is the roundabout trick to introduce time lag by manually\n", +"//shifting the respone by 0.7 min\n", +"plot(t,y)\n", +"xtitle('Ex-6.3(Fig 6.9)','Time(min)','y(t)/KM');\n", +"a=legend('Nonlinear regression','Smiths method','FOPTD',position=4);\n", +"a.font_size=2;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;\n", +"//====NON-LINEAR REGRESSION====//\n", +"//Now that we have the response data and also taking the word from the book that\n", +"//simulation from Excel/Matlab is identical to experimental data, we can actually \n", +"//take this simulation and use it for showing regression\n", +"//So our experimental data is y(1)\n", +"//For nonlinear regression we define a cost function which we have to minimize\n", +"function err=experiment(tau)\n", +" s=%s;\n", +" G=syslin('c',1/((tau(1)*s+1)*(tau(2)*s+1)));\n", +" t=0:0.1:20;\n", +" response=csim('step',t,G);\n", +" err=sum((response-y(1,:)).^2);\n", +"endfunction\n", +"//f is value of cost function, g is gradient of cost function, \n", +"//ind is just a dummy variable required by optim function\n", +"function [f,g,ind]=cost(tau,ind) \n", +" f=experiment(tau)\n", +" g=numdiff(experiment,tau)\n", +"endfunction\n", +"x0=[3 1]'; //Initial guess for optimization routine\n", +"[cost_opt, tau_opt]=optim(cost,x0)\n", +"mprintf('\n Optimization using optim function done successfully \n')\n", +"mprintf('\n From nonlinear regression \n tau1=%f min\n tau2=%f min \n',[tau_opt]')\n", +"//Formatted output...\n", +"mprintf('\n tau1(min) tau2(min) Sum of squares')\n", +"mprintf('\n Smith %f %f %f',3.81,0.84,0.0769)\n", +"mprintf('\n First Order\n(delay=0.7min) %f %s %f',4.60,'--',0.0323)\n", +"mprintf('\n Excel and Matlab %f %f %f \n',3.34,1.86,0.0057)" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 6.5: Estimation_of_model_parameters.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 6.5\n", +"disp('Example 6.5')\n", +"k=0:10';\n", +"yk=[0 0.058 0.217 0.360 0.488 0.600 0.692 0.772 0.833 0.888 0.925]';\n", +"Y=yk(2:$);\n", +"n=length(Y);\n", +"yk_1=[yk(1:$-1)];\n", +"yk_2=[0;yk(1:$-2)];\n", +"uk_1=ones(n,1);\n", +"uk_2=[0;uk_1(1:$-1)];\n", +"X=[yk_1 yk_2 uk_1 uk_2];\n", +"Bhat=inv(X'*X)*X'*Y;//Equation 6-10\n", +"//a1, a2, b1, b2 are components of Bhat for linear regression\n", +"K_lin=(Bhat(3)+Bhat(4))/(1-Bhat(1)-Bhat(2)); //Equation 6-42\n", +"//====NON-LINEAR REGRESSION====//\n", +"//Now that we have the response data we can do non-linear regression through\n", +"//the transfer function approach. Total no. of variables to be regressed are\n", +"//three: K, tau1, tau2.\n", +"//For nonlinear regression we define a cost function which we have to minimize\n", +"function err=experiment(tau)\n", +" K=tau(3);tau1=tau(1);tau2=tau(2);\n", +" t=k';\n", +" response=tau(3)*(1-(tau1*exp(-t/tau1)-tau2*exp(-t/tau2))/(tau1-tau2))\n", +" err=sum((response-[yk]).^2);\n", +"endfunction\n", +"//f is value of cost function, g is gradient of cost function, \n", +"//ind is just a dummy variable required by optim function\n", +"function [f,g,ind]=cost(tau,ind) \n", +" f=experiment(tau)\n", +" g=numdiff(experiment,tau)\n", +"endfunction\n", +"x0=[1 3 1]'; //Initial guess for optimization routine\n", +"[cost_opt, tau_opt]=optim(cost,x0)\n", +"mprintf('\n Optimization using optim function done successfully \n')\n", +"mprintf('\n From nonlinear regression \n tau1=%f \n tau2=%f min\n K=%f min \n',[tau_opt]')\n", +"//Now we have to get discrete ARX model parameters from transfer function parameters \n", +"//For this we use Equation nos.: 6-32 to 6-36\n", +"deltat=1;taua=0;tau1=tau_opt(1);tau2=tau_opt(2);K=tau_opt(3);\n", +"a1=exp(-deltat/tau1)+exp(-deltat/tau2);\n", +"a2=-exp(-deltat/tau1)*exp(-deltat/tau2);\n", +"b1=K*(1+(taua-tau1)/(tau1-tau2)*exp(-deltat/tau1)-(taua-tau2)/(tau1-tau2)*exp(-deltat/tau2));\n", +"b2=K*(exp(-deltat*(1/tau1+1/tau2))+(taua-tau1)/(tau1-tau2)*exp(-deltat/tau2)-(taua-tau2)/(tau1-tau2)*exp(-deltat/tau1));\n", +"mprintf('\n Linear Regression Non-Linear Regression')\n", +"mprintf('\n %s %f %20f',['a1';'a2';'b1';'b2';'K'],[[Bhat;K_lin] [a1;a2;b1;b2;K]])\n", +"yL_hat=X*Bhat;\n", +"yN_hat=X*[a1;a2;b1;b2];\n", +"mprintf('\n \n n y yL_hat yN_hat')\n", +"mprintf('\n %f %f %f %f',[1:10]',yk(2:$),yL_hat,yN_hat)\n", +"mprintf('\n \n Note that values of coefficients for non-linear regression \n are different...\n", +"from that of linear regression, but the\n output is the same \n...\n", +"hence this shows that the coefficients need not be unique....\n", +"\n the coefficients for nonlinear regression given in book and from this scilab code\n...\n", +" both are correct')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 6.6: Step_test_of_distillation_column.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 6.6\n", +"disp('Example 6.6')\n", +"mprintf('\n It is not possible to fit Model 1 or \n plot it because experimental data...\n", +" has not been given in the book. \n Hence we simply plot Model 2,3,4 \n')\n", +" \n", +"//Model 2\n", +"a=[3.317 -4.033 2.108 0.392 ]'\n", +"b=[-0.00922 0.0322 -0.0370 0.0141]';\n", +"uk=[ones(120,1)]; //Inputs-step at t=1 min\n", +"yk=zeros(120,1); //Outputs initialization\n", +"for k=5:120\n", +" yk(k)=a(1)*yk(k-1)+a(2)*yk(k-2)+a(3)*yk(k-3)+a(4)*yk(k-4)...\n", +" +b(1)*uk(k-1)+b(2)*uk(k-2)+b(3)*uk(k-3)+b(4)*uk(k-4);\n", +"end\n", +"//Model 2 trial with transfer function\n", +"//a=-flipdim([-1 3.317 -4.033 2.108 0.392 ]',1);\n", +"//b=flipdim([-0.00922 0.0322 -0.0370 0.0141]',1);\n", +"//\n", +"//Gz=poly(b,'z','coeff')/poly(a,'z','coeff');\n", +"//u=ones(120,1);\n", +"//Gz=syslin('d',Gz);\n", +"//yk=flts(u',Gz)\n", +"//Although the code is correct, the values given in the book for the coeffs\n", +"//of the ARX model are wrong and hence the system diverges(blows up)\n", +"mprintf('Although the code is correct, the values \n given in the book for the coeffs of model 2...\n", +"\n of the ARX model are wrong and hence the system diverges(blows up)')\n", +"//Model 3\n", +"s=%s;\n", +"G3=0.082/(7.95*s+1);//We have to add delay of 44.8 min\n", +"//Model 4\n", +"G4=0.088*(1-12.2*s)/(109.2*s^2+23.1*s+1);//We have to add delay of 25.7 min\n", +"G=syslin('c',[G3;G4]);\n", +"t=[0:0.1:120]';\n", +"y=csim('step',t,G);\n", +"y(1,:)=[zeros(1,448) y(1,1:$-448)]\n", +"y(2,:)=[zeros(1,257) y(2,1:$-257)]\n", +"plot2d(t,y')\n", +"xtitle('Ex-6.6','Time(min)','y(t)');\n", +"a=legend('Model-3','Model-4',position=4);\n", +"a.font_size=2;\n", +"a=get('current_axes');b=a.title;b.font_size=5;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/8-Control_System_Instrumentation.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/8-Control_System_Instrumentation.ipynb new file mode 100644 index 0000000..1aebe35 --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/8-Control_System_Instrumentation.ipynb @@ -0,0 +1,127 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 8: Control System Instrumentation" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 8.2: Pump_and_heat_exchanger_system.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 8.2\n", +"disp('Example 8.2')\n", +"//Eqn 8-6\n", +"//Pump characteristics\n", +"q=0:0.1:240;\n", +"Phe=30*(q/200).^2;\n", +"plot2d(q,Phe,rect=[0,0,240,40]);\n", +"xgrid()\n", +"xtitle('Fig 8.13 Pump characteristics','q,gal/min','P,psi')\n", +"scf();\n", +"q=200;//Flow rate in gal/min\n", +"Phe=30*(q/200).^2;\n", +"Pv=40-Phe; //Eqn 8-8\n", +"//(a)\n", +"l=0.5;Pv=10;\n", +"Cv=q/l/sqrt(Pv); \n", +"mprintf('(a) The value of coefficient Cv is %f',Cv)\n", +"//plotting valve characteristic curve\n", +"l=[0:0.01:0.8]';\n", +"n=length(l);\n", +"Cv=125;\n", +"function y=valve_1(q)\n", +" Pv=40-30*(q/200).^2;\n", +"y=Cv*l.*sqrt(Pv)-q;\n", +"endfunction\n", +"[q_valve1,f1]=fsolve(200*ones(n,1),valve_1); //200*ones(n,1) is the initial guess for q\n", +"plot2d(l,q_valve1);\n", +"//(b)\n", +"q=200*110/100; //110% flow rate\n", +"Phe=30*(q/200).^2;\n", +"Pv=40-Phe; //Eqn 8-8\n", +"l=1;\n", +"Cv=q/sqrt(Pv)/l;\n", +"mprintf('\n(b) The value of coefficient Cv is %f',Cv)\n", +"//We use Cv=115;\n", +"Cv=115;\n", +"l=[0.2:0.01:0.9]';\n", +"n=length(l);\n", +"R=50;\n", +"function y=valve_2(q)\n", +" Pv=40-30*(q/200).^2;\n", +" y=[R^(l-1)]*Cv.*sqrt(Pv)-q;\n", +"endfunction\n", +"[q_valve2,f2]=fsolve(150*ones(n,1),valve_2);\n", +"plot2d(l,q_valve2,style=2)\n", +"//(c)\n", +"Cv=1.2*115;\n", +"mprintf('\n(c) The value of coefficient Cv is %f',Cv)\n", +"l=[0.2:0.01:0.9]';\n", +"n=length(l);\n", +"R=50;\n", +"function y=valve_3(q)\n", +" Pv=40-30*(q/200).^2;\n", +" y=[R^(l-1)]*Cv.*sqrt(Pv)-q;\n", +"endfunction\n", +"[q_valve3,f3]=fsolve(linspace(60,200,n)',valve_3); //Initial guess has to be smart for each valve, \n", +"//since we want near linear profile we can give a linear initial guess\n", +"plot2d(l,q_valve3,style=3)\n", +"//(d)\n", +"Cv=0.8*115;\n", +"mprintf('\n(c) The value of coefficient Cv is %f',Cv)\n", +"l=[0.2:0.01:0.9]';\n", +"n=length(l);\n", +"R=50;\n", +"function y=valve_4(q)\n", +" Pv=40-30*(q/200).^2;\n", +" y=[R^(l-1)]*Cv.*sqrt(Pv)-q;\n", +"endfunction\n", +"[q_valve4,f4]=fsolve(linspace(60,200,n)',valve_4); //Initial guess has to be smart for each valve, \n", +"//since we want near linear profile we can give a linear initial guess\n", +"plot2d(l,q_valve4,style=4,rect=[0,0,1,240])\n", +"xtitle('Ex-8.2 Installed valve characteristics','$l$','q gal/min');\n", +"a=legend('Valve 1, linear Cv=125','Valve 2, Equal% Cv=115','Valve 3, Equal% Cv=138','Valve 4, Equal% Cv=92',position=4);\n", +"a.font_size=2;\n", +"a=get('current_axes');b=a.title;b.font_size=3;c=a.x_label;c.font_size=5;\n", +"c=a.y_label;c.font_size=5;" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Process_Dynamics_And_Controls_by_D_E_Seborg/9-Process_Safety_and_Process_Control.ipynb b/Process_Dynamics_And_Controls_by_D_E_Seborg/9-Process_Safety_and_Process_Control.ipynb new file mode 100644 index 0000000..68169b4 --- /dev/null +++ b/Process_Dynamics_And_Controls_by_D_E_Seborg/9-Process_Safety_and_Process_Control.ipynb @@ -0,0 +1,128 @@ +{ +"cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 9: Process Safety and Process Control" + ] + }, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 9.1: Liquid_surge_system.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 9.1\n", +"disp('Example 9.1')\n", +"q1=12;//cub ft/min\n", +"q2=6;\n", +"q3=13;\n", +"A=%pi*3^2;//ft^2\n", +"delta_t=10;//min\n", +"delta_h_max=delta_t*(q1+q2-q3)/A;\n", +"mprintf('Alarm should be at least %f ft below top of tank',delta_h_max)" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 9.2: Abnormal_event_in_distillation_column.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 9.2\n", +"disp('Example 9.2')\n", +"mu=[0.5 0.8 0.2]; //population means of z y x\n", +"S=[0.01 0.020 0.005];//population std dev of z y x\n", +"z=[0.485]; //steady state values\n", +"y=[0.825];\n", +"x=0.205;\n", +"F=4;D=2;B=2; //flow rates\n", +"Ec=F*z-D*y-B*x;\n", +"disp(Ec,'Ec=')\n", +"sigma_Ec=sqrt(F^2*S(1)^2+D^2*S(2)^2+B^2*S(3)^2)\n", +"disp(sigma_Ec,'sigma_Ec')\n", +"Z=(Ec-0)/sigma_Ec;\n", +"disp(Z,'Z=');\n", +"[P,Q]=cdfnor('PQ',0.120,0,sigma_Ec);\n", +"//Since P is close to 1, we use Q\n", +"Probability=1-2*Q;\n", +"disp(Probability,'Probability of abnormal event=')" + ] + } +, +{ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example 9.3: Reliability_of_flow_control_loop.sce" + ] + }, + { +"cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], +"source": [ +"clear\n", +"clc\n", +"//Example 9.3\n", +"disp('Example 9.3')\n", +"mu=[1.73 0.05 0.49 0.60 0.44]';//failures/yr\n", +"R=exp(-mu);\n", +"P=1-R;\n", +"R_overall=prod(R);\n", +"P_overall=1-R_overall;\n", +"mu_overall=-log(R_overall);\n", +"MTBF=1/mu_overall;\n", +"mprintf('MTBF= %f yr',MTBF)" + ] + } +], +"metadata": { + "kernelspec": { + "display_name": "Scilab", + "language": "scilab", + "name": "scilab" + }, + "language_info": { + "file_extension": ".sce", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "scilab", + "version": "0.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} |