initial commit / add all books

8 files changed, 413 insertions, 0 deletions

diff --git a/1202/CH22/EX22.1/22_1.png b/1202/CH22/EX22.1/22_1.png
new file mode 100755
index 000000000..672cf0288
--- /dev/null
+++ b/1202/CH22/EX22.1/22_1.png
diff --git a/1202/CH22/EX22.1/22_1.sce b/1202/CH22/EX22.1/22_1.sce
new file mode 100755
index 000000000..431391e03
--- /dev/null
+++ b/1202/CH22/EX22.1/22_1.sce
@@ -0,0 +1,53 @@
+clear
+clc
+
+//Example 22.1
+disp('Example 22.1')
+
+//Parameters
+Yxs=0.4;B=0.2;Pm=50;Ki=22;
+a=2.2;mu_m=0.48;Km=1.2;Sf=20;
+
+
+//ODE model
+function ydot=model(t,y,D)
+    X=y(1);S=y(2);P=y(3);
+    
+    Xdot=-D*X+mu(S,P)*X;
+    Sdot=D*(Sf-S)-1/Yxs*mu(S,P)*X;
+    Pdot=-D*P+[a*mu(S,P)+B]*X
+    
+    ydot=[Xdot Sdot Pdot]';
+endfunction
+
+//Rate law
+function mu=mu(S,P)
+    mu=mu_m*(1-P/Pm)*S/(Km+S+S^2/Ki);    
+endfunction
+
+t=0:0.1:100;t0=0;
+y0=[6 5 19.14]';//Initial stable condition
+
+D=0.202*1.1;//10% increase
+y_up = ode(y0, t0, t, list(model,D))
+D=0.202*0.9;//10% decrease
+y_down = ode(y0, t0, t, list(model,D))
+
+subplot(2,1,1);
+plot(t,y_up(1,:),color='red');
+plot(t,y_down(1,:));
+xtitle("$D=0.202\ h^{-1}$","Time(h)","Biomass (g/L)")
+legend("Dilution +10%","Dilution -10%")
+
+subplot(2,1,2);
+t=0:0.1:100;t0=0;
+y0=[6 5 44.05]';//Initial stable condition
+D=0.0389*1.1;//10% increase
+y_up = ode(y0, t0, t, list(model,D))
+D=0.0389*0.9;//10% decrease
+y_down = ode(y0, t0, t, list(model,D))
+
+plot(t,y_up(1,:),color='red');
+plot(t,y_down(1,:))
+xtitle("$D=0.0389\ h^{-1}$","Time(h)","Biomass (g/L)");
+legend("Dilution +10%","Dilution -10%")
diff --git a/1202/CH22/EX22.2/22_2.png b/1202/CH22/EX22.2/22_2.png
new file mode 100755
index 000000000..f7e9edc8c
--- /dev/null
+++ b/1202/CH22/EX22.2/22_2.png
diff --git a/1202/CH22/EX22.2/22_2.sce b/1202/CH22/EX22.2/22_2.sce
new file mode 100755
index 000000000..42c9073e3
--- /dev/null
+++ b/1202/CH22/EX22.2/22_2.sce
@@ -0,0 +1,215 @@
+clear
+clc
+
+//Example 22.2
+disp('Example 22.2')
+//Author: Dhruv Gupta....Aug 4, 2013
+//<dgupta6@wisc.edu>
+
+K=[0.2 0.58 0.35;0.25 1.10 1.3;0.3 0.7 1.2];
+tau=[2 2 2;3 3 3;4 4 4];
+s=%s;
+
+G=K./(1+tau*s);
+
+RGA=K.*inv(K');
+disp(RGA,"RGA=")
+
+//IMC based tuning
+tauC=5;
+Kc=diag(tau/tauC./K);
+mprintf("\n\nThe tauI given in book are wrong\n...
+refer to Table 11.1 for correct formula\n\n")
+tauI=diag(tau)+1; 
+mprintf('\nWe still however use the ones given in book\n');
+
+
+disp(Kc,"Kc=")
+disp(tauI,"tauI=")
+//Refer to Eqns 15-23 and 15-24
+Gc=Kc.*(1+(1)./tauI/s);
+//For the sake of brevity we write Gstar as G
+//We will account for delays in the for loop that we will write
+//Refer to Figure 15.9 Page 295 for details of Smith Predictor
+
+
+//====Making step response models of the continuos transfer functions====//
+Ts=0.1;//Sampling time ie delta_T
+delay=3/Ts;
+N=150/Ts;//Model Order
+s=%s;
+G=syslin('c',diag(matrix(G,1,9)));//Transfer function
+t=0:Ts:N*Ts;
+u_sim=ones(9,length(t));
+//u_sim(:,1:4)=zeros(9,4); //input delay to account for 3 min delay in G
+S=csim(u_sim,t,G)';//generating step response model for real plant
+//plot(t,S);
+S(1,:)=[];
+//Now we have these step response models for each of the transfer functions
+//[S1 S4 S7
+//S2 S5 S8
+//S3 S6 S9]
+
+T=150+delay;//Simulation Run Time in minutes(we add delay because our for loop runs till n-delay)
+n=T/Ts*2+1; //no. of discrete points in our domain of analysis
+//Input initialization T is the Time for simulation
+
+//========Set point as 10=============//
+//p is the controller output
+p=zeros(n,3);
+delta_p=zeros(n,3);
+ytilde=zeros(n,3); //Prediction of Smith Fig 15.9
+e=zeros(n,3); //corrections
+edash=zeros(n,3);
+delta_edash=zeros(n,3);
+ysp=zeros(n,3);
+ysp((n-1)/2+1:n,1)=10*ones(n-((n-1)/2+1)+1,1);
+
+t=-(n-1)/2*Ts:Ts:(n-1)/2*Ts;
+y=zeros(n,3);
+
+for k=(n-1)/2+1:n-delay
+    
+    //Error  e
+    e(k,:)=ysp(k-1,:)-y(k-1,:);
+    
+    //Error edash
+    edash(k,:)=e(k-1,:)-ytilde(k-1,:)+ytilde(k-1-delay,:);
+    //Edash=E-(Y1-Y2)...where Y2 is delayed Y1
+    delta_edash(k,:)=edash(k,:)-edash(k-1,:);
+
+    //Controller calculation----Digital PID----Eqn 7-28 Pg 136 (Velocity form)
+    p(k,:)=p(k-1,:)+[delta_edash(k,:)+edash(k,:)*diag(Ts./tauI)]*diag(Kc);
+    
+    //Limits on manipulated variables
+    p(k,:)=min([(345-180)*ones(1,3);p(k,:)],'r');
+    p(k,:)=max([(105-180)*ones(1,3);p(k,:)],'r');
+    
+    delta_p(k,:)=p(k,:)-p(k-1,:);
+    
+    
+    //Prediction
+    ytilde(k,1)=[S(1:N-1,1);S(1:N-1,4);S(1:N-1,7)]'...
+        *[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)]...
+        +[S(N,1) S(N,4) S(N,7)]*[p(k-N,1);p(k-N,2);p(k-N,3)];
+    ytilde(k,2)=[S(1:N-1,2);S(1:N-1,5);S(1:N-1,8)]'...
+        *[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)]...
+        +[S(N,2) S(N,5) S(N,8)]*[p(k-N,1);p(k-N,2);p(k-N,3)];
+    ytilde(k,3)=[S(1:N-1,3);S(1:N-1,6);S(1:N-1,9)]'...
+        *[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)]...
+        +[S(N,3) S(N,6) S(N,9)]*[p(k-N,1);p(k-N,2);p(k-N,3)];
+    
+    //Output
+    y(k+delay,1)=[S(1:N-1,1);S(1:N-1,4);S(1:N-1,7)]'...
+        *[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)]...
+        +[S(N,1) S(N,4) S(N,7)]*[p(k-N,1);p(k-N,2);p(k-N,3)];
+    y(k+delay,2)=[S(1:N-1,2);S(1:N-1,5);S(1:N-1,8)]'...
+        *[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)]...
+        +[S(N,2) S(N,5) S(N,8)]*[p(k-N,1);p(k-N,2);p(k-N,3)];
+    y(k+delay,3)=[S(1:N-1,3);S(1:N-1,6);S(1:N-1,9)]'...
+        *[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)]...
+        +[S(N,3) S(N,6) S(N,9)]*[p(k-N,1);p(k-N,2);p(k-N,3)];
+end
+
+
+subplot(2,2,1);
+plot(t',y(:,1),'--',t',y(:,2),':',t',y(:,3),'-.',t',ysp(:,1),'-');
+set(gca(),"data_bounds",[0 150 -4 12]); //putting bounds on display
+l=legend("$y1$","$y2$","$y3$",position=1);
+l.font_size=5;
+xtitle("","Time(min)","$y$");
+a=get("current_axes");
+c=a.y_label;c.font_size=5;
+
+
+subplot(2,2,2);
+plot(t',p(:,1),'--',t',p(:,2),':',t',p(:,3),'-.');
+set(gca(),"data_bounds",[-1 150 -40 100]); //putting bounds on display
+l=legend("$p1$","$p2$","$p3$",position=1);
+l.font_size=5;
+xtitle("","Time(min)","$p$");
+a=get("current_axes");
+c=a.y_label;c.font_size=5;
+
+mprintf("Note that there is no overshoot around time=25 mins \n...
+which is in contrast to what is shown in book")
+
+
+//========Now for set point as 50=============//
+
+//p is the controller output
+p=zeros(n,3);
+delta_p=zeros(n,3);
+ytilde=zeros(n,3); //Prediction of Smith Fig 15.9
+e=zeros(n,3); //corrections
+edash=zeros(n,3);
+delta_edash=zeros(n,3);
+ysp=zeros(n,3);
+ysp((n-1)/2+1:n,1)=50*ones(n-((n-1)/2+1)+1,1);
+
+t=-(n-1)/2*Ts:Ts:(n-1)/2*Ts;
+y=zeros(n,3);
+
+for k=(n-1)/2+1:n-delay
+    
+    //Error  e
+    e(k,:)=ysp(k-1,:)-y(k-1,:);
+    
+    //Error edash
+    edash(k,:)=e(k-1,:)-ytilde(k-1,:)+ytilde(k-1-delay,:);
+    //Edash=E-(Y1-Y2)...where Y2 is delayed Y1
+    delta_edash(k,:)=edash(k,:)-edash(k-1,:);
+
+    //Controller calculation----Digital PID----Eqn 7-28 Pg 136 (Velocity form)
+    p(k,:)=p(k-1,:)+[delta_edash(k,:)+edash(k,:)*diag(Ts./tauI)]*diag(Kc);
+    
+    //Limits on manipulated variables
+    p(k,:)=min([(345-180)*ones(1,3);p(k,:)],'r');
+    p(k,:)=max([(105-180)*ones(1,3);p(k,:)],'r');
+    
+    delta_p(k,:)=p(k,:)-p(k-1,:);
+    
+    
+    //Prediction
+    ytilde(k,1)=[S(1:N-1,1);S(1:N-1,4);S(1:N-1,7)]'...
+        *[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)]...
+        +[S(N,1) S(N,4) S(N,7)]*[p(k-N,1);p(k-N,2);p(k-N,3)];
+    ytilde(k,2)=[S(1:N-1,2);S(1:N-1,5);S(1:N-1,8)]'...
+        *[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)]...
+        +[S(N,2) S(N,5) S(N,8)]*[p(k-N,1);p(k-N,2);p(k-N,3)];
+    ytilde(k,3)=[S(1:N-1,3);S(1:N-1,6);S(1:N-1,9)]'...
+        *[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)]...
+        +[S(N,3) S(N,6) S(N,9)]*[p(k-N,1);p(k-N,2);p(k-N,3)];
+    
+    //Output
+    y(k+delay,1)=[S(1:N-1,1);S(1:N-1,4);S(1:N-1,7)]'...
+        *[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)]...
+        +[S(N,1) S(N,4) S(N,7)]*[p(k-N,1);p(k-N,2);p(k-N,3)];
+    y(k+delay,2)=[S(1:N-1,2);S(1:N-1,5);S(1:N-1,8)]'...
+        *[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)]...
+        +[S(N,2) S(N,5) S(N,8)]*[p(k-N,1);p(k-N,2);p(k-N,3)];
+    y(k+delay,3)=[S(1:N-1,3);S(1:N-1,6);S(1:N-1,9)]'...
+        *[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)]...
+        +[S(N,3) S(N,6) S(N,9)]*[p(k-N,1);p(k-N,2);p(k-N,3)];
+end
+
+
+subplot(2,2,3);
+plot(t',y(:,1),'--',t',y(:,2),':',t',y(:,3),'-.',t',ysp(:,1),'-');
+set(gca(),"data_bounds",[0 150 -10 60]); //putting bounds on display
+l=legend("$y1$","$y2$","$y3$",position=1);
+l.font_size=5;
+xtitle("","Time(min)","$y$");
+a=get("current_axes");
+c=a.y_label;c.font_size=5;
+
+
+subplot(2,2,4);
+plot(t',p(:,1),'--',t',p(:,2),':',t',p(:,3),'-.');
+set(gca(),"data_bounds",[-1 150 -100 200]); //putting bounds on display
+l=legend("$p1$","$p2$","$p3$",position=1);
+l.font_size=5;
+xtitle("","Time(min)","$p$");
+a=get("current_axes");
+c=a.y_label;c.font_size=5;
+
diff --git a/1202/CH22/EX22.3/22_3.png b/1202/CH22/EX22.3/22_3.png
new file mode 100755
index 000000000..19e02fb9d
--- /dev/null
+++ b/1202/CH22/EX22.3/22_3.png
diff --git a/1202/CH22/EX22.3/22_3.sce b/1202/CH22/EX22.3/22_3.sce
new file mode 100755
index 000000000..d93f2d307
--- /dev/null
+++ b/1202/CH22/EX22.3/22_3.sce
@@ -0,0 +1,52 @@
+
+clear
+clc
+
+//Example 22.3
+disp('Example 22.3')
+
+//Parameters
+p1=0.028735;p2=0.028344;p3=5.035E-5;V1=12;n=0.0926;
+Ib=15;//basal
+Gb=81;
+
+//Diet function
+function D=D(t)
+    D=9*exp(-0.05*t);
+endfunction
+
+
+//ODE model
+function ydot=model(t,y,U)
+    G=y(1);X=y(2);I=y(3);
+    Gdot=-p1*G-X*(G+Gb)+D(t);
+    Xdot=-p2*X+p3*I;
+    Idot=-n*(I+Ib)+U/V1;
+    ydot=[Gdot Xdot Idot]';
+endfunction
+
+
+t=0:0.1:400;t0=0;
+y0=[0 0 0]';//G,X,I are deviation variables
+
+U=0;
+y = Gb+ode(y0, t0, t, list(model,U))
+subplot(3,1,1);
+plot(t,y(1,:));
+xtitle("","Time(min)","Glucose (mg/L)")
+legend("$U=0\ mU/min$")
+
+U=15;
+y =Gb+ ode(y0, t0, t, list(model,U))
+subplot(3,1,2);
+plot(t,y(1,:));
+xtitle("","Time(min)","Glucose (mg/L)")
+legend("$U=15\ mU/min$")
+
+U=25;
+y = Gb+ode(y0, t0, t, list(model,U))
+subplot(3,1,3);
+plot(t,y(1,:));
+xtitle("","Time(min)","Glucose (mg/L)")
+legend("$U=25\ mU/min$")
+
diff --git a/1202/CH22/EX22.4/22_4.png b/1202/CH22/EX22.4/22_4.png
new file mode 100755
index 000000000..29a133013
--- /dev/null
+++ b/1202/CH22/EX22.4/22_4.png
diff --git a/1202/CH22/EX22.4/22_4.sce b/1202/CH22/EX22.4/22_4.sce
new file mode 100755
index 000000000..4faf109f9
--- /dev/null
+++ b/1202/CH22/EX22.4/22_4.sce
@@ -0,0 +1,93 @@
+clear
+clc
+
+//Example 22.4
+disp('Example 22.4')
+K=[-6 3;12 5];
+tau=[0.67 2;0.67 5];
+s=%s;
+G=K./(1+tau*s);
+delay75=0.75;
+delay1=1;
+RGA=K.*inv(K');
+disp(RGA,"RGA=")
+
+//IMC based tuning
+tauC=[tau(1,1) tau(2,2)];
+Kc=diag(tau./(repmat(tauC,2,1)+[delay75 delay1;delay75 delay1])./K);
+tauI=diag(tau);
+disp(Kc,"Kc=")
+disp(tauI,"tauI=")
+
+//====Making step response models of the continuos transfer functions====//
+Ts=0.05;//Sampling time ie delta_T
+delay75=0.75/Ts;
+delay1=1/Ts;
+N=30/Ts;//Model Order
+s=%s;
+G=syslin('c',diag(matrix(G,1,4)));//Transfer function
+t=0:Ts:N*Ts;
+u_sim=ones(4,length(t));
+u_sim([1 2],1:(delay75))=zeros(2,delay75); //input delay to account for delay in SNP
+u_sim([3 4],1:(delay1))=zeros(2,delay1); //input delay to account for delay in DPM
+S=csim(u_sim,t,G)';//generating step response model for real plant
+//plot(t,S);
+S(1,:)=[];
+//Now we have these step response models for each of the transfer functions
+//[S1 S3
+//S2 S4 
+
+
+
+
+T=50;//Simulation Run Time in minutes
+n=T/Ts*2+1; //no. of discrete points in our domain of analysis
+
+
+//========Set point as -10=============//
+//p is the controller output
+p=zeros(n,2);
+delta_p=zeros(n,2);
+e=zeros(n,2); //errors=(ysp-y) on which PI acts
+ysp=zeros(n,2);
+ysp((n-1)/2+1:n,1)=-10*ones(n-((n-1)/2+1)+1,1);
+
+t=-(n-1)/2*Ts:Ts:(n-1)/2*Ts;
+y=zeros(n,2);
+
+
+for k=(n-1)/2+1:n
+    
+    //Error  e
+    e(k,:)=ysp(k-1,:)-y(k-1,:);
+    delta_e(k,:)=e(k,:)-e(k-1,:);
+
+    //Controller calculation----Digital PID----Eqn 7-28 Pg 136 (Velocity form)
+//    p(k,:)=p(k-1,:)+flipdim([delta_e(k,:)+e(k,:)*diag(Ts./tauI)]*diag(Kc),2);
+    p(k,:)=p(k-1,:)+([delta_e(k,:)+e(k,:)*diag(Ts./tauI)]*diag(Kc));
+    //1-1/2-2 pairing
+    
+    delta_p(k,:)=p(k,:)-p(k-1,:);
+       
+    //Output
+    y(k,1)=[S(1:N-1,1);S(1:N-1,3)]'...
+        *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1)]...
+        +[S(N,1) S(N,3)]*[p(k-N,1);p(k-N,2)];
+    y(k,2)=[S(1:N-1,2);S(1:N-1,4)]'...
+        *[flipdim(delta_p(k-N+1:k-1,1),1);flipdim(delta_p(k-N+1:k-1,2),1);]...
+        +[S(N,2) S(N,4)]*[p(k-N,1);p(k-N,2)];
+end
+
+
+plot(t',y(:,1),'--',t',y(:,2),':',t',ysp(:,1),'-');
+set(gca(),"data_bounds",[0 50 -20 25]); //putting bounds on display
+l=legend("Mean arterial pressure","Cardiac output",position=1);
+//l.font_size=5;
+xtitle("","Time(min)","$y$");
+a=get("current_axes");
+c=a.y_label;c.font_size=5;
+
+mprintf('\nThere is more interaction in the variables\n...
+than the book claims, hence a mismatch between the result\n...
+and the book\n')
+