Code covered by the BSD License

# Marine Automatics

### Mateusz Blonski (view profile)

28 Oct 2010 (Updated )

Simulation models of marine automatics elements library for MATLAB/Simulink

Demo2cs2.m
```% Script for solving kd coefficient under the I2 index for Demo no. 2
% (C) Copyright by Mateusz Blonski, Szczecin 2010

clear;
clc;

% Define variables
syms s kd kdt kp J L R b ku km b2 b1 b0 a3 a2 a1 a0 B3 B2 B1;

% Set the transfer function of controller
Gr=kp;
display('Gr=');
pretty(kp);

% Set the transfer function of motor
Gm=collect(km/(L*J*s^3+(R*J+L*b)*s^2+(R*b+ku*km)*s));
display('Gm=');
pretty(Gm);

% Set the transfer function of internal loop of the system
Gi=collect((Gr*Gm)/(1+Gr*Gm*kd*s));
display('Gi=');
pretty(Gi);

% Set the transfer function of error signal
e=collect((1/(1+Gi*kdt))/s);
display('e=');
pretty(e);

% Calculate the W(s) polynomial
W=collect((b2*s^2+b1*s+b0)*(b2*s^2-b1*s+b0));
display('W=');
pretty(W);

% Define coefficients of I2 index
b2=J*L;
b1=L*b+J*R;
b0=R*b+km*(ku+kd*kp);
a3=b2;
a2=b1;
a1=b0;
a0=kdt*km*kp;
B3=simple(b2^2);
B2=simple(2*b0*b2-b1^2);
B1=simple(b0^2);

% Calculate I2 index
I2=simple((a3*a2*B1-a3*a0*B2+a1*a0*B3)/(2*a3*a0*(a2*a1-a3*a0)));
display('I2=');
pretty(I2);

% Differentiate the I2 index over kd coefficient
eq=simple(diff(I2,kd));
display('eq=');
pretty(eq);

% Solve equation eq with kd as the variable
solv=solve(eq,kd);
display('kd=');
pretty(simple(solv));```