#*-------------------------------------------------------------------
* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
*
* This LIBRARY is free software; you can distribute it and/or modify
* it under the therms of the ALSOC FREE LICENSE as available at
* http://www.enq.ufrgs.br/alsoc.
*
* EMSO is distributed under the terms of the ALSOC LICENSE as
* available at http://www.enq.ufrgs.br/alsoc.
*-----------------------------------------------------------------------
* Author: Jonathan Ospino P.
* $Id: Tuning.mso  2012$
*---------------------------------------------------------------------*#

#* Ideal PID Tuning tool  

Brief description 
------------------

This tool allows the user to obtain the controller tuning parameters according to some of the
most common PID tuning rules found in the literature.In order to be used this block, one must supply
the following parameters:

CASE 1: The parameters of the different elements of the control system are known:
* Process characterization: Gain(K), Time Constant(Tau), and Dead Time (t0; if present). 
* Sensor-Transmitter: Gain(Km) and Time Constant(Taum).
* Valve: Gain(Kf) and Time Constant(Tauf).

CASE 2: The parameters of the FODPDT approximation of the response curve are known:
In this case the process characterization are taken as the parameters of the FODPDT approximation.
For making the algorithm constistent, the parameters of the other elements of control system must be set to ZERO!!! 


What results are reported to the user?
--------------------------------------

According to the type of controller and the parameters specified, the block does the respective 
calculations and  shows the following results for each of the different tuning rules:

* Controller tuning parameters:
	- Controller proportional gain (Kc)
	- Controller integral time (TauI)
	- Controller derivative time (TauD)
	- Controller proportional  band (PB)
	- Controller reset rate (TauI_R)

References
*Corripio and Smith. Principles of Automatic Process Control. 2005
*Luyben  and Luyben. Essentials of Process Control.1997
*Marlin T. Process control : designing processes and control systems for dynamic performance. 2000 

*#

using "types";

Model PID_Tuning

ATTRIBUTES
Pallete=true;
Icon="icon/Tuning";
Info="== Ideal PID Tuning ==

	It computes the tuning parameters of an Ideal PID according to different tuning rules. 
	This tool can be used in the following two cases: 
	
	(1) If the parameters of the basic elements of the control system are known 
	(Gain, time constant, and time delay). 
	
	(2) If the parameters of a FODPDT approximation of the response curve are known.
	 
	In the  first case all the parameters must be specified normally. 
	In the second case, however, the parameters of all the basic elements of the control system, 
	but the process, must be specified as ZERO. In this case the process parameters are taken 
	as the parameters of the FODPDT approximation of the response curve.";


PARAMETERS

# QUALITATIVE PARAMETERS
Controller_Type as Switcher(Valid=["P","PI","PID"],Default="PID");
Lopez_et_al_Criterion as Switcher(Valid=["IAE","ISE","ITAE"],Default="IAE");
Rovira_et_al_Criterion as Switcher(Valid=["IAE","ITAE"],Default="IAE");

#PROCESS
Kp as Real(Brief="<<Process Steady-State Gain>> or <<Gain of the response curve as an approximation to a FODPDT>>",Default=1);
Taup as Real(Brief="<<Process Time Constant>> or <<Time Constant of the response curve as an approximation to a FODPDT>>",Default=1);
t0p as Real(Brief="<<Process Time Delay>> or <<Time Delay of the response curve as an approximation to a FODPDT>>",Default=0);

#SENSOR-TRANSMITTER
Km as Real(Brief="Sensor-Transmitter Gain (NOTE:Set it to ZERO in case of using a Response curve characterization)",Default=1);
Taum as Real(Brief="Sensor-Transmitter Time Constant (NOTE:Set it to ZERO in case of using a Response curve characterization)",Default=0);

#VALVE
Kf as Real(Brief="Valve Steady-State Gain (NOTE:Set it to ZERO in case of using a Response curve characterization)",Default=1);
Tauf as Real(Brief="Valve Time Constant (NOTE:Set it to ZERO in case of using a Response curve characterization)",Default=1);


VARIABLES

# PARAMETERS FOR THE Kcu CALCULATION PROCEDURE USING A FIRST-ORDER PADÉ APPROXIMATION

#* Coefficients of the polynomial used for computing the ultimate parameters
   in case of having the FIRST CASE listed below with the presence of time delay.*#
A as Real(Hidden=true);
B as Real(Hidden=true);
C as Real(Hidden=true);
#* Variables used to store the value of the calculated roots of the wu polynomial*#
wu12 as Real(Hidden=true);
wu22 as Real(Hidden=true);

# PARAMETERS OF THE FODPDT APPROXIMATION
K as Real(Brief="Gain of the FODPDT approximation",Protected=true);
Tau as Real(Brief="Time Constant of the FODPDT approximation",Protected=true);
t0 as Real(Brief="Time Delay of the FODPDT approximation",Protected=true);

# PARAMETERS OF THE LOPEZ'S TUNING RULE 
aL1 as Real(Hidden=true);
bL1 as Real(Hidden=true);
aL2 as Real(Hidden=true);
bL2 as Real(Hidden=true);
aL3 as Real(Hidden=true);
bL3 as Real(Hidden=true);

# PARAMETERS OF THE ROVIRAS'S TUNING RULE 
aR1 as Real(Hidden=true);
bR1 as Real(Hidden=true);
aR2 as Real(Hidden=true);
bR2 as Real(Hidden=true);
aR3 as Real(Hidden=true);
bR3 as Real(Hidden=true);

# Results for each set of Tuning Rules  
ZieglerNichols_ClosedLoop as TunedController(Protected=true);
TyreusLuyben_ClosedLoop as TunedController(Protected=true);
ZieglerNichols_OpenLoop as TunedController(Protected=true);
CohenCoon_OpenLoop as TunedController(Protected=true);
Lopez_et_al_Performance as TunedController(Protected=true);
Rovira_et_al_Performance as TunedController(Protected=true);
Ciancone_Regulatory as TunedController(Protected=true);
Ciancone_Servo as TunedController(Protected=true);

# Parameters used in the stability criteria
Kcu as Real(Brief="Controller Ultimate Gain",Protected=true);
Tu as Real(Brief="Controller Ultimate Period",Protected=true);
wu as Real(Brief="Controller Ultimate frequency",Protected=true);


EQUATIONS

# ********************************************************
# *** 1. Previous calculations to the tuning procedure ***
# ********************************************************

# ********************************************
# *** 1.1. Ultimate parameters calculation *** 
# ********************************************

if abs(Kp*Kf*Km)>0 then 
	## FIRST CASE: All the parameters of the control system are known
	if t0p<=0 then 
		### Systems without time delay 		
		if Taum>0 then 
			#### Considering the Sensor-Transmitter time constant (3 FOD)
			
			"Direct Substitution Method on a 3FOD system"
			wu=sqrt((Taup+Taum+Tauf)/(Taup*Taum*Tauf));
			Kcu=((Tauf*Taup+Tauf*Taum+Taum*Taup)*(Tauf+Taup+Taum)-Taup*Tauf*Taum)/(Tauf*Taup*Taum*Kf*Kp*Km);
			Tu=2*3.141592/wu;
			
			# Parameters of the most complex case
			A=0;
			B=0;
			C=0;
			wu12=0;
			wu22=0;
		else 
			#### Without considering the Sensor-Transmitter time constant (2 FOD)
			
			"Direct Substitution Method on a 2FOD system"
			wu=0;
			Kcu=-1/(Kf*Kp*Km); #It has to be checked it out!!!
			Tu=0;
			
			# Parameters of the most complex case
			A=0;
			B=0;
			C=0;
			wu12=0;
			wu22=0;
		end
	else 
		### Systems with time delay
		if Taum>0 then 
			#### Considering the Sensor-Transmitter time constant (3 FOD + Time Delay)
			
			#*	
			By the Direct substitution method coupled with a First-Order Padé approximation
			we can obtain a Fourth-Degree polynomial of wu of the form:
	
			A*wu^4 + B*wu^2 + C= 0
	
			which can be used for computing wu in a simple way, as follows: 
			*#
			
			"Direct Substitution Method on a 2FOD+1FODPDT system"
			A=-(Tauf*Taup*Taum)*t0p^2/2;
			B=t0p*((t0p*(Tauf+Taup+Taum)+2*(Tauf*Taup+Tauf*Taum+Taum*Taup))/2-(t0p*(Tauf*Taup+Tauf*Taum+Taum*Taup)+2*Taup*Taum*Tauf));
			C=-2*(t0p+Taup+Taum+Tauf);
			wu12=(-B+sqrt(B^2-4*A*C))/(2*A);
			wu22=(-B-sqrt(B^2-4*A*C))/(2*A);
			
			if wu12<=0 then
				if wu22>0 then
					wu=sqrt(wu22);
				else
					wu=0;
				end
			else
				wu=sqrt(wu12); # wu12 is greater than wu22
			end
			
			Kcu=-(t0p*(Tauf*Taup*Taum)/(2*Kp*Kf*Km))*wu^4+((t0p*(Taup+Taum+Tauf)+2*(Taup*Taum+Taum*Tauf+Tauf*Taup))/(2*Kp*Kf*Km))*wu^2-1/(Kp*Kf*Km);
			if wu>0 then
				Tu=2*3.141592/wu;
			else
				Tu=0;
			end
		else 
			#### Without considering the Sensor-Transmitter time constant (TO BE CHECKED !!!)
			
			#*	
			By the Direct substitution method coupled with a First-Order Padé approximation
			we can obtain a Quadratic polynomial of wu of the form:
	
			A*wu^2 + B = 0
	
			which can be used for computing wu in a simple way, as follows:
			*#		
			
			"Direct Substitution Method on a 1FOD+1FODPDT system"
			A=-(Tauf*Taup*t0p+t0p*(t0p*Tauf+t0p*Taup+2*Tauf*Taup)/2);
			B=2*(t0p+Tauf+Taup);
			C=0;
			wu12=0;
			wu22=-B/A; # temporarily saves the result
			
			if wu22>0 then
				wu=sqrt(wu22);
			else
				wu=0;
			end
			
			Kcu=((t0p*Tauf+t0p*Taup+2*Tauf*Taup)*wu^2-2)/(2*Kp*Kf*Km);
			
			if wu>0 then
				Tu=2*3.141592/wu;
			else
				Tu=0;
			end			
		end		
	end
else 
	## SECOND CASE: The parameters of a FODPDT response curve are known
	if t0p>0 then
		### System with time delay	

		"Direct Substitution Method on a FODPDT with a First-Order Padé Approximation"
		wu=(2/t0p)*sqrt(t0p/Taup+1);
		Kcu=(1+2*Taup/t0p)/Kp;
		Tu=2*3.141592/wu;

		# Parameters of the most complex case
		A=0;
		B=0;
		C=0;
		wu12=0;
		wu22=0;
	else
		### System without time delay
		
		"Direct Substitution Method on a FOD"
		wu=0;
		Kcu=-1/Kp; 
		Tu=0;
		#* WARNING: 
			This case doesn't have an ultimate value.
			User should be capable of understanding this by oneself.
			This equations are put here just for equaling the number
			of the equations in both side of the IF-ELSE sentence*#
		
		# Parameters of the most complex case
		A=0;
		B=0;
		C=0;
		wu12=0;
		wu22=0;
	end
end


# *****************************************************************
# *** 1.2. Approximation of a 2FOD+1FODPDT to a FODPDT dynamics ***
# *****************************************************************
if (abs(Kf)<=0) and (abs(Km)<=0) then
	# Case 1: An approximation of the response curve to a FODPDT is provided 
	K=Kp;
	Tau=Taup;
	t0=t0p;	
else
	# Case 2: The parameters Kf,Tauf,Kp,Taup,t0p,Km, and Taum are provided
	# The following simple estimation is considered:
	K=Kf*Kp*Km;	
	Tau=max([Taup,Tauf,Taum]);                      
	t0=sum([Taup,Tauf,Taum])-max([Taup,Tauf,Taum]);
	
end


# ******************************************
# *** 2. Application of the Tuning Rules ***
# ******************************************

# ********************************************
# ** 2.1 ZIEGLER-NICHOLS CLOSED-LOOP METHOD **   TO BE CHECKED!!
# ********************************************
switch Controller_Type
	case "P":
	ZieglerNichols_ClosedLoop.Kc=Kcu/2;                            
	ZieglerNichols_ClosedLoop.TauI=0; # It cannot be calculated
	ZieglerNichols_ClosedLoop.TauD=0; # It cannot be calculated
	case "PI":
	ZieglerNichols_ClosedLoop.Kc=Kcu/2.2;                          
	ZieglerNichols_ClosedLoop.TauI=Tu/1.2;                         
	ZieglerNichols_ClosedLoop.TauD=0; # It cannot be calculated
	case "PID":
	ZieglerNichols_ClosedLoop.Kc=1.25*Kcu/1.7;
	ZieglerNichols_ClosedLoop.TauI=5*Tu/8;
	ZieglerNichols_ClosedLoop.TauD=Tu/10;
end


# ******************************************
# ** 2.2 TYREUS-LUYBEN CLOSED-LOOP METHOD **   
# ******************************************
# Taken from Luyben & Luyben. Essentials of process control. 
switch Controller_Type
	case "P":
	# This method doesn't have correlations for a Only-P controller
	TyreusLuyben_ClosedLoop.Kc=0;   # It cannot be calculated
	TyreusLuyben_ClosedLoop.TauI=0; # It cannot be calculated
	TyreusLuyben_ClosedLoop.TauD=0; # It cannot be calculated
	case "PI":
	TyreusLuyben_ClosedLoop.Kc=Kcu/3.2;
	TyreusLuyben_ClosedLoop.TauI=2.2*Tu;
	TyreusLuyben_ClosedLoop.TauD=0;
	case "PID":
	TyreusLuyben_ClosedLoop.Kc=Kcu/2.2;
	TyreusLuyben_ClosedLoop.TauI=2.2*Tu;
	TyreusLuyben_ClosedLoop.TauD=Tu/6.3;
end


# ******************************************
# ** 2.3 ZIEGLER-NICHOLS OPEN-LOOP METHOD **   
# ******************************************
switch Controller_Type
	case "P":
		if t0>0 then
			ZieglerNichols_OpenLoop.Kc=(1/K)*(Tau/t0);
			ZieglerNichols_OpenLoop.TauI=0; # It cannot be calculated
			ZieglerNichols_OpenLoop.TauD=0; # It cannot be calculated
		else
			ZieglerNichols_OpenLoop.Kc=0;
			ZieglerNichols_OpenLoop.TauI=0; # It cannot be calculated
			ZieglerNichols_OpenLoop.TauD=0; # It cannot be calculated
		end 
	case "PI":
		if t0>0 then
			ZieglerNichols_OpenLoop.Kc=(0.9/K)*(Tau/t0);
			ZieglerNichols_OpenLoop.TauI=3.3*t0;
			ZieglerNichols_OpenLoop.TauD=0; # It cannot be calculated
		else
			ZieglerNichols_OpenLoop.Kc=0;
			ZieglerNichols_OpenLoop.TauI=0;
			ZieglerNichols_OpenLoop.TauD=0; # It cannot be calculated
		end
	case "PID":
		if t0>0 then
			ZieglerNichols_OpenLoop.Kc=(1.2/K)*(Tau/t0);
			ZieglerNichols_OpenLoop.TauI=2.0*t0;
			ZieglerNichols_OpenLoop.TauD=0.5*t0;
		else
			ZieglerNichols_OpenLoop.Kc=0;
			ZieglerNichols_OpenLoop.TauI=0;
			ZieglerNichols_OpenLoop.TauD=0;			
		end
end


# *************************************
# ** 2.4 COHEN-COON OPEN-LOOP METHOD **   
# *************************************
switch Controller_Type
	case "P":
		if t0>0 then
			CohenCoon_OpenLoop.Kc=(Tau/(K*t0))*(1+t0/(3*Tau));
			CohenCoon_OpenLoop.TauI=0; # It cannot be calculated
			CohenCoon_OpenLoop.TauD=0; # It cannot be calculated
		else
			CohenCoon_OpenLoop.Kc=0;
			CohenCoon_OpenLoop.TauI=0; # It cannot be calculated
			CohenCoon_OpenLoop.TauD=0; # It cannot be calculated
		end 
	case "PI":
		if t0>0 then
			CohenCoon_OpenLoop.Kc=(Tau/(K*t0))*(0.9+t0/(12*Tau));
			CohenCoon_OpenLoop.TauI=t0*(30+3*t0/Tau)/(9+20*t0/Tau);
			CohenCoon_OpenLoop.TauD=0; # It cannot be calculated
		else
			CohenCoon_OpenLoop.Kc=0;
			CohenCoon_OpenLoop.TauI=0;
			CohenCoon_OpenLoop.TauD=0; # It cannot be calculated
		end
	case "PID":
		if t0>0 then
			CohenCoon_OpenLoop.Kc=(Tau/(K*t0))*(4/3+t0/(4*Tau));
			CohenCoon_OpenLoop.TauI=t0*(32+6*t0/Tau)/(13+8*t0/Tau);
			CohenCoon_OpenLoop.TauD=4*t0/(11+2*t0/Tau);
		else
			CohenCoon_OpenLoop.Kc=0;
			CohenCoon_OpenLoop.TauI=0;
			CohenCoon_OpenLoop.TauD=0;			
		end
end


# *******************************************************************
# ** 2.5 LOPEZ (REGULATORY): TUNING RULES FOR INTEGRAL PERFORMANCE **  
# *******************************************************************
switch Lopez_et_al_Criterion
	case "IAE":
		switch Controller_Type
			case "P":
			aL1=0.902;
			bL1=-0.985;
			aL2=0; # It does not exist
			bL2=0; # It does not exist
			aL3=0; # It does not exist
			bL3=0; # It does not exist
			Lopez_et_al_Performance.Kc=(aL1/K)*(t0/Tau)^(bL1);
			Lopez_et_al_Performance.TauI=0; # It cannot be calculated
			Lopez_et_al_Performance.TauD=0; # It cannot be calculated
			case "PI":
			aL1=0.984;
			bL1=-0.986;
			aL2=0.608; 
			bL2=0.707; 
			aL3=0; # It does not exist
			bL3=0; # It does not exist
			Lopez_et_al_Performance.Kc=(aL1/K)*(t0/Tau)^(bL1);
			Lopez_et_al_Performance.TauI=(Tau/aL2)*(t0/Tau)^(bL2);
			Lopez_et_al_Performance.TauD=0; # It cannot be calculated
			case "PID":
			aL1=1.435;
			bL1=-0.921;
			aL2=0.878; 
			bL2=0.749; 
			aL3=0.482; 
			bL3=1.137; 
			Lopez_et_al_Performance.Kc=(aL1/K)*(t0/Tau)^(bL1);
			Lopez_et_al_Performance.TauI=(Tau/aL2)*(t0/Tau)^(bL2);
			Lopez_et_al_Performance.TauD=aL3*Tau*(t0/Tau)^(bL3);
		end
	case "ISE":
		switch Controller_Type
			case "P":
			aL1=1.411;
			bL1=-0917;
			aL2=0; # It does not exist
			bL2=0; # It does not exist
			aL3=0; # It does not exist
			bL3=0; # It does not exist
			Lopez_et_al_Performance.Kc=(aL1/K)*(t0/Tau)^(bL1);
			Lopez_et_al_Performance.TauI=0; # It cannot be calculated
			Lopez_et_al_Performance.TauD=0; # It cannot be calculated
			case "PI":
			aL1=1.305;
			bL1=-0.959;
			aL2=0.492; 
			bL2=0.739; 
			aL3=0; # It does not exist
			bL3=0; # It does not exist
			Lopez_et_al_Performance.Kc=(aL1/K)*(t0/Tau)^(bL1);
			Lopez_et_al_Performance.TauI=(Tau/aL2)*(t0/Tau)^(bL2);
			Lopez_et_al_Performance.TauD=0; # It cannot be calculated
			case "PID":
			aL1=1.495;
			bL1=-0.945;
			aL2=1.101; 
			bL2=0.771; 
			aL3=0.560; 
			bL3=1.006; 
			Lopez_et_al_Performance.Kc=(aL1/K)*(t0/Tau)^(bL1);
			Lopez_et_al_Performance.TauI=(Tau/aL2)*(t0/Tau)^(bL2);
			Lopez_et_al_Performance.TauD=aL3*Tau*(t0/Tau)^(bL3);
		end
	case "ITAE":
		switch Controller_Type
			case "P":
			aL1=0.490;
			bL1=-1.084;
			aL2=0; # It does not exist
			bL2=0; # It does not exist
			aL3=0; # It does not exist
			bL3=0; # It does not exist
			Lopez_et_al_Performance.Kc=(aL1/K)*(t0/Tau)^(bL1);
			Lopez_et_al_Performance.TauI=0; # It cannot be calculated
			Lopez_et_al_Performance.TauD=0; # It cannot be calculated
			case "PI":
			aL1=0.859;
			bL1=-0.977;
			aL2=0.674; 
			bL2=0.680; 
			aL3=0; # It does not exist
			bL3=0; # It does not exist
			Lopez_et_al_Performance.Kc=(aL1/K)*(t0/Tau)^(bL1);
			Lopez_et_al_Performance.TauI=(Tau/aL2)*(t0/Tau)^(bL2);
			Lopez_et_al_Performance.TauD=0; # It cannot be calculated
			case "PID":
			aL1=1.357;
			bL1=-0.947;
			aL2=0.842; 
			bL2=0.738; 
			aL3=0.381; 
			bL3=0.995; 
			Lopez_et_al_Performance.Kc=(aL1/K)*(t0/Tau)^(bL1);
			Lopez_et_al_Performance.TauI=(Tau/aL2)*(t0/Tau)^(bL2);
			Lopez_et_al_Performance.TauD=aL3*Tau*(t0/Tau)^(bL3);
		end
end


# ***************************************************************
# ** 2.6 ROVIRA (SERVO): TUNING RULES FOR INTEGRAL PERFORMANCE **   
# ***************************************************************
switch Rovira_et_al_Criterion
	case "IAE":
		switch Controller_Type
			case "P":
			aR1=0; # It does not exist
			bR1=0; # It does not exist
			aR2=0; # It does not exist
			bR2=0; # It does not exist
			aR3=0; # It does not exist
			bR3=0; # It does not exist
			Rovira_et_al_Performance.Kc=0; # It cannot be calculated
			Rovira_et_al_Performance.TauI=0; # It cannot be calculated
			Rovira_et_al_Performance.TauD=0; # It cannot be calculated
			case "PI":
			aR1=0.758; 
			bR1=-0.861; 
			aR2=1.02; 
			bR2=-0.323; 
			aR3=0; # It does not exist
			bR3=0; # It does not exist
			Rovira_et_al_Performance.Kc=(aR1/K)*(t0/Tau)^(bR1); 
			Rovira_et_al_Performance.TauI=Tau/(aR2+bR2*(t0/Tau)); 
			Rovira_et_al_Performance.TauD=0; # It cannot be calculated
			case "PID":
			aR1=1.086; 
			bR1=-0.869; 
			aR2=0.740; 
			bR2=-0.130; 
			aR3=0.348; 
			bR3=0.914; 
			Rovira_et_al_Performance.Kc=(aR1/K)*(t0/Tau)^(bR1); 
			Rovira_et_al_Performance.TauI=Tau/(aR2+bR2*(t0/Tau)); 
			Rovira_et_al_Performance.TauD=aR3*Tau*(t0/Tau)^(bR3);
		end
	case "ITAE":
		switch Controller_Type
			case "P":
			aR1=0; # It does not exist
			bR1=0; # It does not exist
			aR2=0; # It does not exist
			bR2=0; # It does not exist
			aR3=0; # It does not exist
			bR3=0; # It does not exist
			Rovira_et_al_Performance.Kc=0; # It cannot be calculated
			Rovira_et_al_Performance.TauI=0; # It cannot be calculated
			Rovira_et_al_Performance.TauD=0; # It cannot be calculated
			case "PI":
			aR1=0.586; 
			bR1=-0.916; 
			aR2=1.03; 
			bR2=-0.165; 
			aR3=0; # It does not exist
			bR3=0; # It does not exist
			Rovira_et_al_Performance.Kc=(aR1/K)*(t0/Tau)^(bR1); 
			Rovira_et_al_Performance.TauI=Tau/(aR2+bR2*(t0/Tau)); 
			Rovira_et_al_Performance.TauD=0; # It cannot be calculated
			case "PID":
			aR1=0.965; 
			bR1=-0.855; 
			aR2=0.796; 
			bR2=-0.147; 
			aR3=0.308; 
			bR3=0.9292; 
			Rovira_et_al_Performance.Kc=(aR1/K)*(t0/Tau)^(bR1); 
			Rovira_et_al_Performance.TauI=Tau/(aR2+bR2*(t0/Tau)); 
			Rovira_et_al_Performance.TauD=aR3*Tau*(t0/Tau)^(bR3);
		end
end


# ************************************************************
# ** 2.7 CIANCONE (REGULATORY): TUNING RULES BY CORRELATION ** 
# ************************************************************
switch Controller_Type
	case "P":
	Ciancone_Regulatory.Kc=0;
	Ciancone_Regulatory.TauI=0;
	Ciancone_Regulatory.TauD=0;
	case "PI":
	if (t0/(t0+Tau)) <= 1.0 then
		if (t0/(t0+Tau)) >= 0.9 then
			(Ciancone_Regulatory.Kc)*K=((0.3-0.35)/(1.0-0.9))*((t0/(t0+Tau))-0.9)+0.35;
			(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.5-0.55)/(1.0-0.9))*((t0/(t0+Tau))-0.9)+0.55;
			(Ciancone_Regulatory.TauD)/(t0+Tau)=0;
		else
			if (t0/(t0+Tau)) >= 0.8 then
				(Ciancone_Regulatory.Kc)*K=((0.35-0.4)/(0.9-0.8))*((t0/(t0+Tau))-0.8)+0.4;
				(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.55-0.58)/(0.9-0.8))*((t0/(t0+Tau))-0.8)+0.58;
				(Ciancone_Regulatory.TauD)/(t0+Tau)=0;
			else
				if (t0/(t0+Tau)) >= 0.7 then
					(Ciancone_Regulatory.Kc)*K=((0.4-0.5)/(0.8-0.7))*((t0/(t0+Tau))-0.7)+0.5;
					(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.58-0.605)/(0.8-0.7))*((t0/(t0+Tau))-0.7)+0.605;
					(Ciancone_Regulatory.TauD)/(t0+Tau)=0;
				else
					if (t0/(t0+Tau)) >= 0.6 then
						(Ciancone_Regulatory.Kc)*K=((0.5-0.6)/(0.7-0.6))*((t0/(t0+Tau))-0.6)+0.6;
						(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.605-0.65)/(0.7-0.6))*((t0/(t0+Tau))-0.6)+0.65;
						(Ciancone_Regulatory.TauD)/(t0+Tau)=0;
					else
						if (t0/(t0+Tau)) >= 0.5 then
							(Ciancone_Regulatory.Kc)*K=((0.6-0.9)/(0.6-0.5))*((t0/(t0+Tau))-0.5)+0.9;
							(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.65-0.7)/(0.6-0.5))*((t0/(t0+Tau))-0.5)+0.7;
							(Ciancone_Regulatory.TauD)/(t0+Tau)=0;		
						else
							if (t0/(t0+Tau)) >= 0.4 then
								(Ciancone_Regulatory.Kc)*K=((0.9-1.0)/(0.5-0.4))*((t0/(t0+Tau))-0.4)+1.0;
								(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.7-0.75)/(0.5-0.4))*((t0/(t0+Tau))-0.4)+0.75;
								(Ciancone_Regulatory.TauD)/(t0+Tau)=0;
							else
								if (t0/(t0+Tau)) >= 0.3 then
									(Ciancone_Regulatory.Kc)*K=((1.0-1.13)/(0.4-0.3))*((t0/(t0+Tau))-0.3)+1.13;
									(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.75-0.8)/(0.4-0.3))*((t0/(t0+Tau))-0.3)+0.8;
									(Ciancone_Regulatory.TauD)/(t0+Tau)=0;
								else
									if (t0/(t0+Tau)) >= 0.2 then
										(Ciancone_Regulatory.Kc)*K=((1.13-1.15)/(0.3-0.2))*((t0/(t0+Tau))-0.2)+1.15;
										(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.8-0.52)/(0.3-0.2))*((t0/(t0+Tau))-0.2)+0.52;
										(Ciancone_Regulatory.TauD)/(t0+Tau)=0;
									else
										if (t0/(t0+Tau)) >= 0.1 then
											(Ciancone_Regulatory.Kc)*K=((1.15-1.1)/(0.2-0.1))*((t0/(t0+Tau))-0.1)+1.1;
											(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.52-0.25)/(0.2-0.1))*((t0/(t0+Tau))-0.1)+0.25;
											(Ciancone_Regulatory.TauD)/(t0+Tau)=0;
										else
											(Ciancone_Regulatory.Kc)*K=1.1;
											(Ciancone_Regulatory.TauI)/(t0+Tau)=0.25;
											(Ciancone_Regulatory.TauD)/(t0+Tau)=0;
										end
									end
								end
							end
						end
					end
				end
			end
		end
	else
		Ciancone_Regulatory.Kc=0;
		Ciancone_Regulatory.TauI=0;
		Ciancone_Regulatory.TauD=0;
	end
	case "PID":
	if (t0/(t0+Tau)) <= 1.0 then
		if (t0/(t0+Tau)) >= 0.9 then
			(Ciancone_Regulatory.Kc)*K=((0.3-0.35)/(1.0-0.9))*((t0/(t0+Tau))-0.9)+0.35;
			(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.42-0.46)/(1.0-0.9))*((t0/(t0+Tau))-0.9)+0.46;
			(Ciancone_Regulatory.TauD)/(t0+Tau)=((0.05-0.13)/(1.0-0.9))*((t0/(t0+Tau))-0.9)+0.13;
		else
			if (t0/(t0+Tau)) >= 0.8 then
				(Ciancone_Regulatory.Kc)*K=((0.35-0.45)/(0.9-0.8))*((t0/(t0+Tau))-0.8)+0.45;
				(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.46-0.5)/(0.9-0.8))*((t0/(t0+Tau))-0.8)+0.5;
				(Ciancone_Regulatory.TauD)/(t0+Tau)=((0.13-0.2)/(0.9-0.8))*((t0/(t0+Tau))-0.8)+0.2;
			else
				if (t0/(t0+Tau)) >= 0.7 then
					(Ciancone_Regulatory.Kc)*K=((0.45-0.6)/(0.8-0.7))*((t0/(t0+Tau))-0.7)+0.6;
					(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.5-0.55)/(0.8-0.7))*((t0/(t0+Tau))-0.7)+0.55;
					(Ciancone_Regulatory.TauD)/(t0+Tau)=((0.2-0.15)/(0.8-0.7))*((t0/(t0+Tau))-0.7)+0.15;
				else
					if (t0/(t0+Tau)) >= 0.6 then
						(Ciancone_Regulatory.Kc)*K=((0.6-0.7)/(0.7-0.6))*((t0/(t0+Tau))-0.6)+0.7;
						(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.55-0.6)/(0.7-0.6))*((t0/(t0+Tau))-0.6)+0.6;
						(Ciancone_Regulatory.TauD)/(t0+Tau)=((0.15-0.1)/(0.7-0.6))*((t0/(t0+Tau))-0.6)+0.1;
					else
						if (t0/(t0+Tau)) >= 0.5 then
							(Ciancone_Regulatory.Kc)*K=((0.7-0.9)/(0.6-0.5))*((t0/(t0+Tau))-0.5)+0.9;
							(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.6-0.65)/(0.6-0.5))*((t0/(t0+Tau))-0.5)+0.65;
							(Ciancone_Regulatory.TauD)/(t0+Tau)=((0.1-0.075)/(0.6-0.5))*((t0/(t0+Tau))-0.5)+0.075;
						else
							if (t0/(t0+Tau)) >= 0.4 then
								(Ciancone_Regulatory.Kc)*K=((0.9-1.05)/(0.5-0.4))*((t0/(t0+Tau))-0.4)+1.05;
								(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.65-0.675)/(0.5-0.4))*((t0/(t0+Tau))-0.4)+0.675;
								(Ciancone_Regulatory.TauD)/(t0+Tau)=((0.075-0.05)/(0.5-0.4))*((t0/(t0+Tau))-0.4)+0.05;
							else
								if (t0/(t0+Tau)) >= 0.3 then
									(Ciancone_Regulatory.Kc)*K=((1.05-1.1)/(0.4-0.3))*((t0/(t0+Tau))-0.3)+1.1;
									(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.675-0.7)/(0.4-0.3))*((t0/(t0+Tau))-0.3)+0.7;
									(Ciancone_Regulatory.TauD)/(t0+Tau)=((0.05-0.02)/(0.4-0.3))*((t0/(t0+Tau))-0.3)+0.02;
								else
									if (t0/(t0+Tau)) >= 0.2 then
										(Ciancone_Regulatory.Kc)*K=((1.1-1.2)/(0.3-0.2))*((t0/(t0+Tau))-0.2)+1.2;
										(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.7-0.52)/(0.3-0.2))*((t0/(t0+Tau))-0.2)+0.52;
										(Ciancone_Regulatory.TauD)/(t0+Tau)=((0.02-0.0)/(0.3-0.2))*((t0/(t0+Tau))-0.2)+0;
									else
										if (t0/(t0+Tau)) >= 0.1 then
											(Ciancone_Regulatory.Kc)*K=((1.2-1.1)/(0.2-0.1))*((t0/(t0+Tau))-0.1)+1.1;
											(Ciancone_Regulatory.TauI)/(t0+Tau)=((0.52-0.25)/(0.2-0.1))*((t0/(t0+Tau))-0.1)+0.25;
											(Ciancone_Regulatory.TauD)/(t0+Tau)=0;
										else
											(Ciancone_Regulatory.Kc)*K=1.1;
											(Ciancone_Regulatory.TauI)/(t0+Tau)=0.25;
											(Ciancone_Regulatory.TauD)/(t0+Tau)=0;
										end
									end
								end
							end
						end
					end
				end
			end
		end
	else
		Ciancone_Regulatory.Kc=0;
		Ciancone_Regulatory.TauI=0;
		Ciancone_Regulatory.TauD=0;
	end
end


# *******************************************************
# ** 2.8 CIANCONE (SERVO): TUNING RULES BY CORRELATION **   
# *******************************************************
switch Controller_Type
	case "P":
	Ciancone_Servo.Kc=0;
	Ciancone_Servo.TauI=0;
	Ciancone_Servo.TauD=0;
	case "PI":
	if (t0/(t0+Tau)) <= 1.0 then
		if (t0/(t0+Tau)) >= 0.9 then
			(Ciancone_Servo.Kc)*K=((0.63-0.65)/(1.0-0.9))*((t0/(t0+Tau))-0.9)+0.65;
			(Ciancone_Servo.TauI)/(t0+Tau)=((0.43-0.45)/(1.0-0.9))*((t0/(t0+Tau))-0.9)+0.45;
			(Ciancone_Servo.TauD)/(t0+Tau)=0;
		else
			if (t0/(t0+Tau)) >= 0.8 then
				(Ciancone_Servo.Kc)*K=0.65;
				(Ciancone_Servo.TauI)/(t0+Tau)=((0.45-0.5)/(0.9-0.8))*((t0/(t0+Tau))-0.8)+0.5;
				(Ciancone_Servo.TauD)/(t0+Tau)=0;
			else
				if (t0/(t0+Tau)) >= 0.7 then
					(Ciancone_Servo.Kc)*K=0.65;
					(Ciancone_Servo.TauI)/(t0+Tau)=((0.5-0.57)/(0.8-0.7))*((t0/(t0+Tau))-0.7)+0.57;
					(Ciancone_Servo.TauD)/(t0+Tau)=0;
				else
					if (t0/(t0+Tau)) >= 0.6 then
						(Ciancone_Servo.Kc)*K=((0.65-0.7)/(0.7-0.6))*((t0/(t0+Tau))-0.6)+0.7;
						(Ciancone_Servo.TauI)/(t0+Tau)=((0.57-0.64)/(0.7-0.6))*((t0/(t0+Tau))-0.6)+0.64;
						(Ciancone_Servo.TauD)/(t0+Tau)=0;
					else
						if (t0/(t0+Tau)) >= 0.5 then
							(Ciancone_Servo.Kc)*K=((0.7-0.85)/(0.6-0.5))*((t0/(t0+Tau))-0.5)+0.85;
							(Ciancone_Servo.TauI)/(t0+Tau)=((0.64-0.7)/(0.6-0.5))*((t0/(t0+Tau))-0.5)+0.7;
							(Ciancone_Servo.TauD)/(t0+Tau)=0;		
						else
							if (t0/(t0+Tau)) >= 0.4 then
								(Ciancone_Servo.Kc)*K=((0.85-0.9)/(0.5-0.4))*((t0/(t0+Tau))-0.4)+0.9;
								(Ciancone_Servo.TauI)/(t0+Tau)=((0.7-0.78)/(0.5-0.4))*((t0/(t0+Tau))-0.4)+0.78;
								(Ciancone_Servo.TauD)/(t0+Tau)=0;
							else
								if (t0/(t0+Tau)) >= 0.3 then
									(Ciancone_Servo.Kc)*K=((0.9-1.0)/(0.4-0.3))*((t0/(t0+Tau))-0.3)+1.0;
									(Ciancone_Servo.TauI)/(t0+Tau)=((0.78-0.86)/(0.4-0.3))*((t0/(t0+Tau))-0.3)+0.86;
									(Ciancone_Servo.TauD)/(t0+Tau)=0;
								else
									if (t0/(t0+Tau)) >= 0.2 then
										(Ciancone_Servo.Kc)*K=((1.0-1.05)/(0.3-0.2))*((t0/(t0+Tau))-0.2)+1.05;
										(Ciancone_Servo.TauI)/(t0+Tau)=((0.86-0.95)/(0.3-0.2))*((t0/(t0+Tau))-0.2)+0.95;
										(Ciancone_Servo.TauD)/(t0+Tau)=0;
									else
										if (t0/(t0+Tau)) >= 0.1 then
											(Ciancone_Servo.Kc)*K=((1.05-1.1)/(0.2-0.1))*((t0/(t0+Tau))-0.1)+1.1;
											(Ciancone_Servo.TauI)/(t0+Tau)=((0.95-0.75)/(0.2-0.1))*((t0/(t0+Tau))-0.1)+0.55;
											(Ciancone_Servo.TauD)/(t0+Tau)=0;
										else
											(Ciancone_Servo.Kc)*K=1.1;
											(Ciancone_Servo.TauI)/(t0+Tau)=0.75;
											(Ciancone_Servo.TauD)/(t0+Tau)=0;
										end
									end
								end
							end
						end
					end
				end
			end
		end
	else
		Ciancone_Servo.Kc=0;
		Ciancone_Servo.TauI=0;
		Ciancone_Servo.TauD=0;
	end
	case "PID":
	if (t0/(t0+Tau)) <= 1.0 then
		if (t0/(t0+Tau)) >= 0.9 then
			(Ciancone_Servo.Kc)*K=0.36;
			(Ciancone_Servo.TauI)/(t0+Tau)=((0.42-0.46)/(1.0-0.9))*((t0/(t0+Tau))-0.9)+0.46;
			(Ciancone_Servo.TauD)/(t0+Tau)=((0.04-0.14)/(1.0-0.9))*((t0/(t0+Tau))-0.9)+0.14;
		else
			if (t0/(t0+Tau)) >= 0.8 then
				(Ciancone_Servo.Kc)*K=((0.36-0.4)/(0.9-0.8))*((t0/(t0+Tau))-0.8)+0.4;
				(Ciancone_Servo.TauI)/(t0+Tau)=((0.46-0.5)/(0.9-0.8))*((t0/(t0+Tau))-0.8)+0.5;
				(Ciancone_Servo.TauD)/(t0+Tau)=((0.14-0.22)/(0.9-0.8))*((t0/(t0+Tau))-0.8)+0.22;
			else
				if (t0/(t0+Tau)) >= 0.7 then
					(Ciancone_Servo.Kc)*K=((0.4-0.5)/(0.8-0.7))*((t0/(t0+Tau))-0.7)+0.5;
					(Ciancone_Servo.TauI)/(t0+Tau)=((0.5-0.55)/(0.8-0.7))*((t0/(t0+Tau))-0.7)+0.55;
					(Ciancone_Servo.TauD)/(t0+Tau)=((0.22-0.16)/(0.8-0.7))*((t0/(t0+Tau))-0.7)+0.16;
				else
					if (t0/(t0+Tau)) >= 0.6 then
						(Ciancone_Servo.Kc)*K=((0.5-0.6)/(0.7-0.6))*((t0/(t0+Tau))-0.6)+0.6;
						(Ciancone_Servo.TauI)/(t0+Tau)=((0.55-0.6)/(0.7-0.6))*((t0/(t0+Tau))-0.6)+0.6;
						(Ciancone_Servo.TauD)/(t0+Tau)=((0.16-0.12)/(0.7-0.6))*((t0/(t0+Tau))-0.6)+0.12;
					else
						if (t0/(t0+Tau)) >= 0.5 then
							(Ciancone_Servo.Kc)*K=((0.6-0.7)/(0.6-0.5))*((t0/(t0+Tau))-0.5)+0.7;
							(Ciancone_Servo.TauI)/(t0+Tau)=((0.6-0.68)/(0.6-0.5))*((t0/(t0+Tau))-0.5)+0.68;
							(Ciancone_Servo.TauD)/(t0+Tau)=((0.12-0.08)/(0.6-0.5))*((t0/(t0+Tau))-0.5)+0.08;
						else
							if (t0/(t0+Tau)) >= 0.4 then
								(Ciancone_Servo.Kc)*K=((0.7-0.9)/(0.5-0.4))*((t0/(t0+Tau))-0.4)+0.9;
								(Ciancone_Servo.TauI)/(t0+Tau)=((0.68-0.8)/(0.5-0.4))*((t0/(t0+Tau))-0.4)+0.8;
								(Ciancone_Servo.TauD)/(t0+Tau)=((0.08-0.06)/(0.5-0.4))*((t0/(t0+Tau))-0.4)+0.06;
							else
								if (t0/(t0+Tau)) >= 0.3 then
									(Ciancone_Servo.Kc)*K=((0.9-1.04)/(0.4-0.3))*((t0/(t0+Tau))-0.3)+1.04;
									(Ciancone_Servo.TauI)/(t0+Tau)=((0.8-0.87)/(0.4-0.3))*((t0/(t0+Tau))-0.3)+0.87;
									(Ciancone_Servo.TauD)/(t0+Tau)=((0.06-0.04)/(0.4-0.3))*((t0/(t0+Tau))-0.3)+0.04;
								else
									if (t0/(t0+Tau)) >= 0.2 then
										(Ciancone_Servo.Kc)*K=((1.04-1.07)/(0.3-0.2))*((t0/(t0+Tau))-0.2)+1.07;
										(Ciancone_Servo.TauI)/(t0+Tau)=((0.87-0.94)/(0.3-0.2))*((t0/(t0+Tau))-0.2)+0.94;
										(Ciancone_Servo.TauD)/(t0+Tau)=((0.04-0.02)/(0.3-0.2))*((t0/(t0+Tau))-0.2)+0.02;
									else
										if (t0/(t0+Tau)) >= 0.1 then
											(Ciancone_Servo.Kc)*K=((1.07-1.1)/(0.2-0.1))*((t0/(t0+Tau))-0.1)+1.1;
											(Ciancone_Servo.TauI)/(t0+Tau)=((0.94-0.74)/(0.2-0.1))*((t0/(t0+Tau))-0.1)+0.74;
											(Ciancone_Servo.TauD)/(t0+Tau)=((0.02-0.0)/(0.2-0.1))*((t0/(t0+Tau))-0.1)+0.0;
										else
											(Ciancone_Servo.Kc)*K=1.1;
											(Ciancone_Servo.TauI)/(t0+Tau)=0.74;
											(Ciancone_Servo.TauD)/(t0+Tau)=0;
										end
									end
								end
							end
						end
					end
				end
			end
		end
	else
		Ciancone_Servo.Kc=0;
		Ciancone_Servo.TauI=0;
		Ciancone_Servo.TauD=0;
	end
end

end

Model TunedController
		
	VARIABLES
	
	Kc as Real(Brief="Controller Gain",Protected=true);
	TauI as Real(Brief="Controller Integral Time",Protected=true);
	TauD as Real(Brief="Controller Derivative Time",Protected=true);
	PB as Real(Brief="Controller Proportional Band",Protected=true);
	TauI_R as Real(Brief="Controller reset rate",Protected=true);


	EQUATIONS
	
	if Kc>0 then
		PB=100/Kc;
	else
		PB=0;
	end
	
	if TauI>0 then
		TauI_R=1/TauI;
	else
		TauI_R=0;
	end
	
end