// Functions       : FitPID_tv
// Author          : Josef Schuetzenberger - Austria
// License         : Public Domain
// Version         : 1.0
// Date            : 2020-09-09
// Release         : http://fileexchange.scilab.org/toolboxes/585000

// Supported Scilab Version: 6.1

// DESCRIPTION
// Numerically computing (fitting) a Transfer Function to a measured response(y) 
// of an unknown system. The input function(x) could be a step, impulse or whatever.
//
// NOTES
// A step includes a wide range of frequencies and is therefore a good choice.
// You may have to subtract the steady state first.
 
// PARAMETERS:
//   Inputs:
//     x: Vector of an impulse or step or slope or...
//     y: Vector of a response to the step-function or an impulse or... (Same length than x)
//     tv: The given derivative time 
//
//   Outputs:   
//     G: Transfer Function
//     kp: is the proportional gain
//     tn: is the integral time
//     kv: is the derivative gain
//     tv: is the given derivative time
//     of: offset
 
clear() //Clear Variables
close() //Close graphic figures, progression bars or wait bars

function [G,kp,tn,kv,tv,of]=FitPID_tv(x,y,tv)//Transfer-function of a PID-T1 Controller (1+s+s^2)/(s+s^2) (tv is given)
    A=[];x1=0;y1=0;x2=0             
    for k=1:length(x), 
        x1=x1+x(k)
        x2=x2+x1
        y1=y1+y(k)
        y(k)= y(k)+y1/tv
        A($+1,:)=[x2,x1,x(k),1]
    end
    X= A \ y'
    of= X(4)  //offset
    tn= X(2)/X(1)-tv 
    kp= X(1)*tv*tn
    kv= X(3)/kp-1
    G= syslin('c',(X(1)+X(2)*%s+X(3)*%s^2)*1,1/tv*%s+%s^2) 
endfunction

// Demo 

// Recorded data from a csv-file
M=csvRead("FitPID_tv_Set.csv", [],[], "double")
y=M(:,3)
xd=M(:,1)-M(:,2)
o=y(1);y=y-o            // subtract steady state
//The function would not be able to find the correct parameters without the pre-set of tv. 
//There is probably not enough information in the recordings. (Sampling time, Quantization)
[G,kp,tn,kv,tv,of]=FitPID_tv(xd',y',60) //preset tv=60
disp(G,"fitted settings of the controller kp="+string(kp)+" tn="+string(tn)+"s tv="+string(tv)+"s kv="+string(kv))
disp("actual settings of the controller kp="+string(6)+" tn="+string(52)+"s tv="+string(60)+"s kv="+string(1))
yr= flts(xd',dscr(G,1)) // calc response from fitted transfer funktion
yr=yr+o                 // add steady state
Mp=[M(:,1),M(:,2),M(:,3),yr']
plot(Mp)
title("Fitted response of a PID-Controller to a sticky valve (tv given)", "fontsize",3)
xlabel("Time (s)")
ylabel("Magnitude in (%)")
legend("X-Input","W-Input (Setpoint)","Y-Output" ,"Fitted Y-Output",4,%t);
xgrid(1)
xstring( 950, 62,msprintf('%s %3.1f ',"kp=",kp))
xstring( 950, 60,msprintf('%s %3.0f %s',"tn=",tn,"s"))
xstring( 950, 58,msprintf('%s %3.0f %s',"tv=",tv,"s"))
xstring( 950, 56,msprintf('%s %3.1f \n',"kv=",kv))
xstring( 950, 50, prettyprint(G))
