TractionDrive_P

function TractionDrive_P
%See: "Different energetic techniques for modelling traction drives"
%Zanasi R., Geitner G.-H., Bouscayrol A., Lhomme W.; Proc. of the 2008
%International Conference ElectrIMACS, Quebec; CD

assignin('base','Cf',0.006);      %Filter capacitance [F=As/V]
Rf=0.01;                          %Filter resistance [Ohm]
assignin('base','Rf',Rf);         %
assignin('base','Lf',0.0009);     %Filter inductance [H=Vs/A]
assignin('base','Lfd1',0.5);      %DC machine 1 field inductance [H=Vs/A]
assignin('base','Lfd2',0.5);      %DC machine 2 field inductance [H=Vs/A]
Rfd1=2;                           %DC machine 1 field resistance [Ohm]
assignin('base','Rfd1',Rfd1);        
assignin('base','Rfd2',2);        %DC machine 2 field resistance [Ohm]
assignin('base','Larm1',0.00054); %DC machine 1 armature inductance [H=Vs/A]
assignin('base','Larm2',0.00054); %DC machine 2 armature inductance [H=Vs/A]
assignin('base','Rarm1',0.025);   %DC machine 1 armature resistance [Ohm]
assignin('base','Rarm2',0.025);   %DC machine 2 armature resistance [Ohm]
assignin('base','Kdcm1',0.077);   %DC machine 1 torque coefficient [Vs/A]
assignin('base','Kdcm2',0.077);   %DC machine 2 torque coefficient [Vs/A]
     %K_gear=8.75; R_wheel=0.483 [m]; m_b1 = mb2 = K_gear/R_wheel;
assignin('base','mb1',18.12);     %bogie ratio 1 [m^-1]
assignin('base','mb2',18.12);     %bogie ratio 2 [m^-1]
assignin('base','M',15000);       %Chassis mass [Kg]
assignin('base','F0',800);        %Resistance force parameter F0 [N]
assignin('base','ar1',0.05);      %Resistance force parameter ar1 [m/s^2]
assignin('base','ar',30);         %Resistance force parameter ar [Ns/m]
assignin('base','br',4);          %Resistance force parameter br [Ns^2/m^2]
assignin('base','alpha',0);       %Resistance force parameter alpha [rad]
assignin('base','gc',9.81);       %Gravity coefficient [m/s^2]
     %Nominal(rated) and maximal values
Udc=750;                          %Input direct voltage [V]
assignin('base','Udc',Udc);       %
I_fd_nom=24;                      %nominal field current [A]
I_arm_nom=440;                    %nominal armature current [A]
     %definition of initial values
ifd10=I_fd_nom;                   %Initial field current DC machine 1 [A]   
assignin('base','ifd10',ifd10);   %
ifd20=I_fd_nom;                   %Initial field current DC machine 2 [A]
assignin('base','ifd20',ifd20);   %
iarm0=0;                          %Initial armature current DC machines [A]
assignin('base','iarm0',iarm0);   %
I_M_0='I_M_0=[ifd10; iarm0; ifd20];';%Initial current vector [A]
evalin('base',I_M_0);             %
i_f_0=(ifd10+iarm0+ifd20)*0+3;    %Initial filter current [A]
assignin('base','i_f_0',i_f_0);   %
u_f_0=Udc-Rf*i_f_0;               %Initial filter voltage [V]
assignin('base','u_f_0',u_f_0);   %
assignin('base','vsub0',0);       %Initial velocity [m/s]                                           %
     %definition of matrices and vectors of the elements
MR0='M_R=';                       %Define resistance matrix
   MR1='[Rfd1      0         0 ;';
   MR2='   0  Rarm1+Rarm2    0 ;';
   MR3='   0       0       Rfd2 ];';
evalin('base',[MR0 MR1 MR2 MR3]); %Berechnung im MATLAB workspace
ML0='M_L=';                       %Define inductance matrix
   ML1='[Lfd1      0         0 ;';
   ML2='   0  Larm1+Larm2    0 ;';
   ML3='   0       0       Lfd2 ];';
evalin('base',[ML0 ML1 ML2 ML3]); %Computation in MATLAB workspace
Mb='M_b=[mb1; mb2];';             %Define bogie vector
evalin('base',Mb);                %Computation in MATLAB workspace
     %Switching computation
m_cp1=(Rfd1*I_fd_nom)/(Udc-Rf*(I_arm_nom+2*I_fd_nom)); %Control computation
%m_cp2=0.3;                                            %
%m_cp3=m_cp1;                                          %
Ts=0.5*0.001;                     %Switching period in sec.
Tc=Ts;
Tfin=6;                           %Simulation period in sec.
assignin('base','Tfin',Tfin);     %
Tcp1=2.3;                         %Computation start
Tcp2=Tcp1+10;
tcp=(0:Tc:Tfin)';
sl_ar=0.23;                       %rising coeff. 1 for mcp2
sl2_ar=0.04;                      %rising coeff. 2 for mcp2
sl_fd=0.009;
Laenge=length(tcp);
mtcp1=zeros(Laenge,1); mtcp2=mtcp1; %preallocation
for ii=1:Laenge
    if tcp(ii)<Tcp1                 %valid till Tcp sec.!
        mtcp1(ii)=m_cp1;            %constant
        mtcp2(ii)=sl_ar*tcp(ii);    %linear rising
    elseif tcp(ii)<Tcp2
        mtcp1(ii)=m_cp1-sl_fd*(tcp(ii)-Tcp1)^0.5; %non-lin decay rate
        if mtcp1(ii)<0.02 mtcp1(ii)=0.02; end     %lower limit
           mtcp2(ii)=sl_ar*Tcp1+sl2_ar*(tcp(ii)-Tcp1); %linear rising
        if mtcp2(ii)>1    mtcp2(ii)=1;    end     %upper limit
    end
end
assignin('base','tcp',tcp);             %
assignin('base','mtcp1',mtcp1);         %
assignin('base','mtcp2',mtcp2);         %
m_cp='m_cp=[tcp mtcp1 mtcp2 mtcp1];';   % Switching functions
evalin('base',m_cp);                    %

txt='"TractionDrive"';
disp(['Parameter for Bond Graph example ' txt ' loaded!']);