기후 변화, 방사능 유출같은 세계적인 문제들로 인해, 그리고 에너지 소모를 줄이기 위해 많은 공학자들이 세로운 제어기 설계 방법을 찾으려고 노력하고 있습니다.
그러한 설계 방법을 그린 엔지니어링 green engineering이라고 부르며, 그린 엔지니어링의 목표는 오염을 줄이고, 사람에게 유해를 줄이도록 하는 제품을 만드는것이 목표라 할수 있겠습니다. 이러한 그린 엔지니어링의 원칙에 따라서 피드벡 제어 시스템의 유용함을 보여드리고자 합니다.
지구 온난화 문제와 오염을 줄이기 위해서는 환경 모니터링 시스템을 늘리고, 품질도 높여야만하는데요. 한 가지의 예시로 외부 환경을 관측하는 무선으로 이동하면서 감지하는 플랫폼들이 있습니다. 보내지는 전력과 전압의 변화 등을 측정하는것도 그러한 예시라고 할수 있겠습니다.
수많은 그린 엔지니어링 시스템들을 다루기 위해서는 신중하게 전류와 전압을 보아야만 하는데요. 전력 변환기는 전력 공급망 상에서 흘러가는 다양한 크기의 전류를 측정하고 모니터링을 할수 있어야 하며, 여기서 사용되는 센서들은 측정치가 제어 시스템이 적절한 동작을 하는데 필요한 정보들을 제공하므로 제어 시스템의 핵심 요소라고 할수 있습니다.
그린 엔지니어링에서의 제어시스템의 역활은 자동화와 정밀화 수준이 높아만 가면서 세계적으로 큰 역활을 하게 될겁니다. 그래서 이후에는 풍력 발전기같은 그린 엔지니어링의 제어나 태양의 밝기가 변화하는 사이에서도 전력 생산량을 최대로 하는 태양광 발전기 피드벡 제어 시스템을 설계해보겠습니다.
제어 시스템으로 개발할수 있는 다른 재밌는 것으로 IoT가 있습니다. 여기서 IoT란 Internet of Things의 약어로 많은 전자 장치들이 연결된것을 말합니다. IoT 네트워크 상에서 수백만의 장치들은 각자의 컴퓨터가 내장되어 인터넷에 연결 됩니다. 이러한 연결된 장비들을 제어해보는것은 제어 공학도에게는 신나는 일이라 할수 있을것 같습니다.
그래서 제어 공학은 재미있고, 도전할만한 분야라고 생각합니다. 제어 시스템은 본래 여러 분야가 혼합된 학문으로, 엔지니어링 커리큘럼 상에서 핵심역활을 한다고 할수 있겠습니다. 이 분야에서는 수학적 배경 지식이 많이 필요하므로, 정리와 증명과 같이 이론적인 부분들도 다루긴 할겁니다. 하지만 최종 목표는 현실에서 동작가능한 제어기를 설계하는 것이므로, ad-hoc(임기 응변)적인 방법과 직관력으로 피드백 제어기를 설계해 보겠습니다.
앞으로 배우는데 있어서 가장 중요한 것은 과거에 다루었던 문제와 이에 대답을 다시 찾아보는것이라 할수 있겠습니다. 그래서 학생들이 수많은 문제들과 해결방안 그리고, 수십년 전에 있어왔던 해결 방안에 대해서 배워보아야 합니다. 옛날의 학습 방법에서 학생은 이러한 문제를 고민하지 않고 단순히 과제들을 풀기만 했었습니다. 하지만 이 과정을 통해서 배우는 사람들이 이론적인 어려움을 줄이고 창의성과 즐거움을 찾길 희망하고 있습니다.
이제 피드벡 제어 이론의 구조에 대해서 살펴보고, 다양한 흥미로운 주제들을 드리고자 합니다. 그래서 여러분들이 피드벡 제어 시스템 이론과 실습에 도움이 되었으면 좋겠습니다.
1. 하나의 전달함수가 주어질때 부분분수 전개 함수 residue로 부분 분수들의 분자와 분모들을 얻습니다.
2. 각 부분분수들을 ilaplace 함수로 역변환 시킨후 변수 after에다가 다 담아줬습니다.
syms s
num = [1 5];
den = [1 6 11 6];
before = tf(num, den)
after = 0;
[r, p, k] = residue(num, den)
% r 분자, p, 분모항
for i = 1:length(r)
F = r(i)/(s - p(i))
after = after + ilaplace(F);
end
after
이번 장에서는 연속 폐루프를 이용해서 횡방향과 종방향 오토파일럿을 설계하였습니다. 횡방향 오토파일럿은 내부 루프로 롤 자세 유지 루프를 외부 루프로 방위각 유지 루프로 이루어집니다.
종방향 오토파일럿은 고도에 따라 더 복잡하게 구성되는데 피치 자세 유지 루프는 모든 영역에서 사용됩니다. 이륙 영역에서는 비행체는 풀 스로틀 상태가 되며 이륙 피치각을 고정하도록 제어하게 됩니다. 상승 영역에서는 비행체는 풀 스로틀로, 피치, 대기 속도 유지 오토 파일럿 루프로 대기속도를 제어합니다. 하강 영역에서는 상승 영역과 비슷한데 대신 비행체에는 최소 스로틀을 줍니다. 고도 유지 영역에서 고도는 피치를 이용한 고도 오토파일럿 루프로 제어되고, 대기속도는 스로틀을 이용한 대기속도 오토파일럿 루프로 제어됩니다.
- 피치 각 제한 $e_{\theta}^{max}$, 댐핑 계수 $\zeta_{\theta}$, $\zeta_{h}$, $\zeta_{V}$, $\zeta_{V2}$와 고유 주파수 $\omega_{n_{v}}$, 고도 루프와 피치 이용 대기속도 루프의 분할 주파수 $W_h$, $W_{V_2}$
계산 고유 주파수들
- 식 (6.21)로 피치 루프 $\omega_{n_{\theta}}$ 고유 주파수 계산. $\omega_{n_{h}}$ = $\omega_{n_{\theta}}$ / $W_h$를 사용하여 고도 루프, $\omega_{n_{V2}}$ = $\omega_{n_{\theta}}$ / $W_{V_2}$를 사용하여 피치 이용 대기속도 루프의 고유 주파수 계산.
% compute trim conditions using 'mavsim_chap5_trim.slx'
% nominal airspeed P.Va0 specified above with aircraft parameters
gamma = 0*pi/180; % desired flight path angle (radians)
R = 10000000000; % desired radius (m) - use (+) for right handed orbit,
% (-) for left handed orbit
Va = 25;
% set initial conditions
x0 = [];
% specify which states to hold equal to the initial conditions
ix = [];
% specify initial inputs
u0 = [...
0;... % 1 - delta_e
0;... % 2 - delta_a
0;... % 3 - delta_r
1;... % 4 - delta_t
];
% specify which inputs to hold constant
iu = [];
% define constant outputs
y0 = [...
Va;... % 1 - Va
0;... % 2 - alpha
0;... % 3 - beta
];
% specify which outputs to hold constant
iy = [1,3];
% define constant derivatives
dx0 =
if R~=Inf, dx0(9) = Va*cos(gamma)/R; end % 9 - psidot
% specify which derivaties to hold constant in trim algorithm
idx = [3; 4; 5; 6; 7; 8; 9; 10; 11; 12];
% compute trim conditions
[x_trim,u_trim,y_trim,dx_trim] = trim('mavsim_trim',x0,u0,y0,ix,iu,iy,dx0,idx);
% check to make sure that the linearization worked (should be small)
norm(dx_trim(3:end)-dx0(3:end))
% P.u_trim = u_trim;
% P.x_trim = x_trim;
% set initial conditions to trim conditions
% initial conditions
MAV.pn0 = 0; % initial North position
MAV.pe0 = 0; % initial East position
MAV.pd0 = -200; % initial Down position (negative altitude)
MAV.u0 = x_trim(4); % initial velocity along body x-axis
MAV.v0 = x_trim(5); % initial velocity along body y-axis
MAV.w0 = x_trim(6); % initial velocity along body z-axis
MAV.phi0 = x_trim(7); % initial roll angle
MAV.theta0 = x_trim(8); % initial pitch angle
MAV.psi0 = x_trim(9); % initial yaw angle
MAV.p0 = x_trim(10); % initial body frame roll rate
MAV.q0 = x_trim(11); % initial body frame pitch rate
MAV.r0 = x_trim(12); % initial body frame yaw rate
trim 시뮬링크 모델이랑 트림 계산 파일이 저렇게 주어졌다.
xdot은 다음과 같으므로
값들은 아래와 같이 대입하면 될거같다
구현 하긴했는데
compute_trim에서 자꾸 dx0이랑 시뮬링크 블록상 연속 상태 크기랑 틀렸다고 에러뜬다.
한참 삽질하다보니 mav_dynamics에서 오일러각대신 쿼터니언을 쓰다보니 상태가 13개가 쓰여서 그렇더라
다시 챕터 3의 mav_dynamics를 오일러각으로 바꿧다.
function [sys,x0,str,ts,simStateCompliance] = mav_dynamics(t,x,u,flag,MAV)
switch flag
%%%%%%%%%%%%%%%%%%
% Initialization %
%%%%%%%%%%%%%%%%%%
case 0
[sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes(MAV);
%%%%%%%%%%%%%%%
% Derivatives %
%%%%%%%%%%%%%%%
case 1
sys=mdlDerivatives(t,x,u,MAV);
%%%%%%%%%%
% Update %
%%%%%%%%%%
case 2
sys=mdlUpdate(t,x,u);
%%%%%%%%%%%
% Outputs %
%%%%%%%%%%%
case 3
sys=mdlOutputs(t,x);
%%%%%%%%%%%%%%%%%%%%%%%
% GetTimeOfNextVarHit %
%%%%%%%%%%%%%%%%%%%%%%%
case 4
sys=mdlGetTimeOfNextVarHit(t,x,u);
%%%%%%%%%%%%%
% Terminate %
%%%%%%%%%%%%%
case 9
sys=mdlTerminate(t,x,u);
%%%%%%%%%%%%%%%%%%%%
% Unexpected flags %
%%%%%%%%%%%%%%%%%%%%
otherwise
DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));
end
% end sfuntmpl
%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes(MAV)
%
% call simsizes for a sizes structure, fill it in and convert it to a
% sizes array.
%
% Note that in this example, the values are hard coded. This is not a
% recommended practice as the characteristics of the block are typically
% defined by the S-function parameters.
%
sizes = simsizes;
sizes.NumContStates = 12;
sizes.NumDiscStates = 0;
sizes.NumOutputs = 12;
sizes.NumInputs = 6;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1; % at least one sample time is needed
sys = simsizes(sizes);
%
% initialize the initial conditions
%
x0 = [...
MAV.pn0;...
MAV.pe0;...
MAV.pd0;...
MAV.u0;...
MAV.v0;...
MAV.w0;...
MAV.phi0;...
MAV.theta0;...
MAV.psi0;...
MAV.p0;...
MAV.q0;...
MAV.r0;...
];
%
% str is always an empty matrix
%
str = [];
%
% initialize the array of sample times
%
ts = [0 0];
% Specify the block simStateCompliance. The allowed values are:
% 'UnknownSimState', < The default setting; warn and assume DefaultSimState
% 'DefaultSimState', < Same sim state as a built-in block
% 'HasNoSimState', < No sim state
% 'DisallowSimState' < Error out when saving or restoring the model sim state
simStateCompliance = 'UnknownSimState';
% end mdlInitializeSizes
%
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%
function sys=mdlDerivatives(t,x,uu, MAV)
pn = x(1);
pe = x(2);
pd = x(3);
u = x(4);
v = x(5);
w = x(6);
phi = x(7);
theta= x(8);
psi = x(9);
p = x(10);
q = x(11);
r = x(12);
fx = uu(1);
fy = uu(2);
fz = uu(3);
ell = uu(4);
m = uu(5);
n = uu(6);
cphi = cos(phi);
ctheta = cos(theta);
cpsi = cos(psi);
sphi = sin(phi);
stheta = sin(theta);
spsi = sin(psi);
pndot = u*ctheta*cpsi -v*(sphi*stheta*cpsi - cphi*spsi) +w*(cphi*stheta*cpsi + sphi*spsi);
pedot = u*ctheta*spsi + v*(sphi*stheta*spsi + cphi*cpsi) + w*(cphi*stheta*spsi - sphi*cpsi);
pddot = -stheta*u + v*sphi*ctheta + w*cphi *ctheta;
udot = r*v - q*w + 1/MAV.mass*fx;
vdot = p*w - r*u + 1/MAV.mass*fy;
wdot = q*u - p*v + 1/MAV.mass*fz;
phidot = p + sphi*tan(theta)*q + r*cphi*tan(theta);
thetadot = q*cphi -sphi*r;
psidot = sphi/ctheta*q + cphi/ctheta*r;
pdot = MAV.Gamma1*p*q-MAV.Gamma2*q*r+MAV.Gamma3*ell+MAV.Gamma4*n;
qdot = MAV.Gamma5*p*r-MAV.Gamma6*(p^2-r^2) +1/MAV.Jy*m;
rdot = MAV.Gamma7*p*q - MAV.Gamma1*q*r + MAV.Gamma4*ell + MAV.Gamma8*n;
sys = [pndot; pedot; pddot; udot; vdot; wdot; phidot; thetadot; psidot; pdot; qdot; rdot];
% end mdlDerivatives
%
%=============================================================================
% mdlUpdate
% Handle discrete state updates, sample time hits, and major time step
% requirements.
%=============================================================================
%
function sys=mdlUpdate(t,x,u)
sys = [];
% end mdlUpdate
%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
%
function sys=mdlOutputs(t,x)
y = [...
x(1);
x(2);
x(3);
x(4);
x(5);
x(6);
x(7);
x(8);
x(9);
x(10);
x(11);
x(12);
];
sys = y;
% end mdlOutputs
%
%=============================================================================
% mdlGetTimeOfNextVarHit
% Return the time of the next hit for this block. Note that the result is
% absolute time. Note that this function is only used when you specify a
% variable discrete-time sample time [-2 0] in the sample time array in
% mdlInitializeSizes.
%=============================================================================
%
function sys=mdlGetTimeOfNextVarHit(t,x,u)
sampleTime = 1; % Example, set the next hit to be one second later.
sys = t + sampleTime;
% end mdlGetTimeOfNextVarHit
%
%=============================================================================
% mdlTerminate
% Perform any end of simulation tasks.
%=============================================================================
%
function sys=mdlTerminate(t,x,u)
sys = [];
% end mdlTerminate
다음 값이 주어질때
트림 값
아래의 트림 시뮬레이션 모델에서
트림 계산 입력이 아래와 같이 주어질때
gamma = 2*pi/180; % desired flight path angle (radians)
R = 5000; % desired radius (m) - use (+) for right handed orbit,
% (-) for left handed orbit
Va = 100;
시뮬링크는 트림 상태를 계산하는 루틴을 제공하는데, help trim 으로 사용방법을 확인할 수 있습니다. 5.3장에서 봤던것 처럼 $V_a^*$, $\gamma^*$, $R^*$가 주어지면, $x^*$, $u^*$를 찾는것이 목표가 됩니다. $\dot{x^*}$ = f($x^*$, $u^*$ 로 여기서 x, u는 식 (5.17)과 (5.18)로 정의 됩니다. $\dot{x^*}$는 식 (5.21)에서 주어지며 f(x, u)는 식 (5.1) - (5.12)로 정의됩니다.
트림 명령의 형식은
[x, u, y, dx] = trim('sys', x0, u0, y0, ix, iu, uy, dx0, idx);
여기서 x는 계산된 트림상태 $x^*$, u는 계산됨 트림 입력 $u^*$, y는 계산된 트림 출력 $y^*$, dx는 계산된 상태의 미분 $\dot{x^*}$, 시뮬링크 모델은 sys.mdl에 명시되는데 모델의 상태들은 서브시스템 내 상태로 되어있고, 입력과 출력은 시뮬링크에서 inports, outports로 정의됩니다.
그림 F.1은 비행체 트림 계산하는 시뮬링크 모델로 입력은 4개의 inports로 서보 명령인 delta_e, delta_a, delta_r, delta_t 가 있습니다. 이 블록의 상태들은 시뮬링크 상태로 우리는 $\zeta$ = ($p_n$, $p_e$, $p_d$, u, u, w, $\phi$, $\theta$, $\psi$, p, q, r)$^T$가 됩니다. 출력은 세개의 outports이며 대기속도 $V_a$, 받음각 $\alpha$, $\beta$이며 이 아웃풋을 트림 명령에다가 전달하여 $V_a$ = $V_a^*$가 유지되도록 하겠습니다.
러더도 사용가능하면 $\beta^*$=0으로 유지시키는 트림 명령을 주어 정상선회 하도록 할수 있씁니다. 러더가 없는 경우 0상태로 만들수 없습니다.
트림 계산에서 문제는 시스템의 비선형 방정식을 풀어야 하는데, 많은 경우 시뮬링크 트림 명령에는 초기화 추정 값으로 상태 x0, 입력 u0, 출력 y0, 상태 미분 dx0이 필요합니다. 우리가 일부 상태나 입력, 출력 또는 미분 상태같은 시작 상태를 알아서 초기 상태를 쓸수 있으면, 사용가능한 값들을 인덱스 벡터 ix, iu, iy, idx로 전달할 수 있습니다.
우리의 상태가 아래와 같다면
이렇게 작성할 수 있습니다.
비슷하게 초기 상태, 입력, 출력들을 다음처럼 정의할수 있습니다.
그림 F.1 트림과 선형 상태 공간 모델 계산을 위한 시뮬링크 다이어그램
F.2 트림 수치 계산 Numerical Computation of Trim
시뮬링크 말고 다른환경에서 시뮬레이션이 개발된다면, 트림 루틴을 따로 작성해야합니다. 이번 장에서는 어떻게 하는지 간단히 설명하고 $V_a^*$, $\gamma^*$, $R^*$ 파라미터들로 상승 선회 트림 제어에 대해 살펴보고, 트림 찾기 알고리즘에서 입력으로 사용하겠습니다. 입력 파라미터 $V_a^*$, $\gamma^*$, $R^*$와 $\alpha$, $\beta$, $\phi$로 트림 상태와 입력 구하는 계산과정을 살펴보겠습니다. 주어진 $V_a^*$, $\gamma^*$, $R^*$값으로 트림값 $\alpha^*$, $\beta^*$, $\phi^*$을 찾아낸다면 이제 트림 상태와 입력같을 구할수 있게됩니다.
첫 단계는 상태 변수들을 $V_a^*$, $\gamma^*$, $R^*$로, 입력을 $\alpha^*$, $\beta^*$, $\phi^*$로 나타내는데 $V_a^*$, $\gamma^*$, $R^*$는 사용자 지정 입력이 되어 알고리즘은 최적의 트림 상태를 구할수 있는 $\alpha^*$, $\beta^*$, $\phi^*$ 계산합니다. 이 값으로 트림 상태 $x^*$, $u^*$를 찾는데 사용됩니다.
동체 좌표계 속도 $u^*$, $v^*$, $w^*$
식 (2.7)에서 동체 좌표계상 속도들은 $V_a^*$, $\alpha^*$, $\beta^*$ 를 이용해서 나타낼 수있습니다.
피치각 $\theta^*$
비행경로각의 정의로 다음과 같이 정리할 수있습니다.
각 변화율 p, q, r
각 변화율은 식 (3.2)로 오일러각으로 나타낼수 있습니다.
엘리베이터 $\delta_e$
$p^*$, $q^*$, $r^*$, $\delta_e^*$에 대해서 식 (5.11)을 풀면
에일러론 $\delta_a$과 러더 $\delta_r$
에일러론과 러더 명령은 식 (5.10), 식(5.12)을 풀어서 구할 수 있습니다.
F.2.1 트림 알고리즘
모든 관심변수와 제어입력은 $V_a^*$, $\gamma^*$, $R^*$,$\alpha^*$, $\beta^*$, $\phi^*$입니다. 트림 알고리즘의 입력은 $V^*$, $\gamma^*$, $R^*$이며, $\alpha^*$, $\beta^*$, $\phi^*$값을 찾기위해 다음의 최적화 알고리즘을 사용합니다.
이는 다음장에서 설명할 경사 하강 알고리즘으로 수치적으로 구할 수있습니다. 트림 알고리즘은 알고리즘 13에서 정리하였습니다.
function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes(MAV)
%
% call simsizes for a sizes structure, fill it in and convert it to a
% sizes array.
%
% Note that in this example, the values are hard coded. This is not a
% recommended practice as the characteristics of the block are typically
% defined by the S-function parameters.
%
sizes = simsizes;
sizes.NumContStates = 13;
sizes.NumDiscStates = 0;
sizes.NumOutputs = 12;
sizes.NumInputs = 6;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1; % at least one sample time is needed
sys = simsizes(sizes);
%
% initialize the initial conditions
%
x0 = [...
MAV.pn0;...
MAV.pe0;...
MAV.pd0;...
MAV.u0;...
MAV.v0;...
MAV.w0;...
MAV.e0;...
MAV.e1;...
MAV.e2;...
MAV.e3;...
MAV.p0;...
MAV.q0;...
MAV.r0;...
];
%
% str is always an empty matrix
%
str = [];
%
% initialize the array of sample times
%
ts = [0 0];
% Specify the block simStateCompliance. The allowed values are:
% 'UnknownSimState', < The default setting; warn and assume DefaultSimState
% 'DefaultSimState', < Same sim state as a built-in block
% 'HasNoSimState', < No sim state
% 'DisallowSimState' < Error out when saving or restoring the model sim state
simStateCompliance = 'UnknownSimState';
% end mdlInitializeSizes
mdlOutputs으로 draw aircraft에서 사용할 변수 12개를 지정해서 출력해주면 된다.
이번장은 매트랩/시뮬링크 환경에 익숙한 사람들을 대상으로 진행하겠습니다. 다세한 정보는 매트랩/시뮬링크 문서를 참조해주시기 바랍니다. 시뮬링크는 상미분 방정식 사이 연결을 푸는대 필수적인 도구입니다.
시뮬링크 각각의 블록들이 다음의 구조대로 되어있다고 가정해봅시다. 여기서 $x_c$ $\in$ $\mathbb{R}$은 초기 조건이 $x_{c0}$인 연속 상태이고, $x_d$ $\in$ $\mathbb{R}^{n_d}$는 초기 조건이 $x_{d0}$인 이산 상태, u $\in$ $\mathbb{R}^m$는 블록 입력, y $\in$ $\mathbb{R}^p$는 블록의 출력이며, t는 시뮬레이션 경과 시간을 의미합니다.
s 함수는 다양한 방법으로 정의 가능한 시뮬링크 도구로 여기서 2가지 방법인 레벨 1 m file s함수와 C파일 s 함수에 대해 간단하게 살펴보겠습니다. c 파일 s 함수는 c코드로 컴파일되서 m 파일 s함수 보다 빠르게 실해오딥니다.
D.1 예시 : 2차 미분 방정식
이 장에서는 표준 2차 전달함수로 시스템을 구현하는 과정을 살펴보겠는데 여기서 레벨 1 m 파일 s 함수와 c 파일 s 함수를 사용하겠습니다.
우선 식 (D.4)를 상태 공간상에 표현하겠습니다. 제어 표준 형태 control canonical form로 나타내면
D.1.1 레벨 1 m 파일 s 함수
식 d.5, d.6을 시스템으로 구현하는 m 파일 s 함수 코드는 아래와 같습니다.
function [sys,x0,str,ts] = second_order_m(t,x,u,flag, zeta,wn)
switch flag,
case 0,
[sys,x0,str,ts]=mdlInitializeSizes;
% initialize block
case 1,
sys=mdlDerivatives(t,x,u,zeta,wn);
% define xdot = f(t,x,u)
case 3,
sys=mdlOutputs(t,x,u,wn);
% define xup = g(t,x,u)
otherwise,
sys = [];
end
%==========================================================%
function [sys,x0,str,ts]=mdlInitializeSizes
sizes = simsizes;
sizes.NumContStates = 2;
sizes.NumDiscStates = 0;
sizes.NumOutputs = 1;
sizes.NumInputs = 1;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;
sys = simsizes(sizes);
x0 = [0; 0]; % define initial conditions
str = []; % str is always an empty matrix
% initialize the array of sample times
ts = [0 0]; % continuous sample time
%==========================================================%
function xdot=mdlDerivatives(t,x,u,zeta,wn)
xdot(1) = −2*zeta*wn*x(1) − wn^2*x(2) + u;
xdot(2) = x(1);
%==========================================================%
function y=mdlOutputs(t,x,u,wn)
y = wn^2*x(2);
여기서 첫 줄은 m 파일 함수를 정의하는 부분이 됩니다. 이 함수의 입력은 항상 경과 시간 t, 상태 x(연속 상태와 이산 상태를 연결한), 입력 u, 와 플래그는 기본적으로 있어야 하며, 사용자 정의 입력 파라미터를 추가할 수도 있습니다. 이번 경우에는 $\zeta$와 $\omega_n$을 사용하겠씁니다.
시뮬링크 엔진은 s 함수를 호출하고 t, x, u, flag를 파라미터로 전달시키는데 flag == 0인경우 시뮬링크 엔진은 sys 구조체를 반환합니다. 여기서 sys 구조체는 초기 상태 x0, 빈 문자열 str, 샘플 타임배열인 ts로 구서오딥니다.
flag == 1인 경우 시뮬링크 엔진은 f(t, x, u) 함수를 반환하고, flag == 2일때는 g(t, x, u), flag==3 일때는 h(t, x, u)를 반환합니다.
스위치 문은 flag를 이용해서 적절한 함수를 호출하게 되는데, 초기화 구문을 살펴보면 연속 상태의 수, 이산 상태의 수, 출력 수, 입력의 수를 정의 합니다. 그리고, 출력을 입력으로 피드백 시킬것인지의 여부를 지정할 수도 있습니다.
