본문으로 바로가기
본 강좌에 사용되는 MATLAB은 버젼 7.9.0 (R2009b)을 대상으로 합니다.

동역학 방정식 !

이떤 사물(혹은 시스템)을 역학적으로 해석해서 미분방정식의 형태로 표현한 것을 동역학 방정식이라고 합니다. 실제로는 참 복잡한 형태를 가지게 되지만, 우리는 그저 간단히 어떻게 동역학을 유도하는지는 생략하교(^^) 어떻게 MATLAB에서 시뮬레이션할 것인지를 이야기해보겠습니다. 이미 이전에 카트형 진자시스템에 대한 이야기를 했었는데요.

이때 사용한 동역학방정식

을 이용하겠습니다. 일단 위 식을 그대로 simulink로 꾸미면, algebraric loop가 있다는 warnning을 보게 됩니다. 왜냐면, 자기자신을 참조(ddot_x, ddot_theta)하기 때문이지요. 그래서 최고차 미분항으로 정리할 필요가 있습니다.

간단히 행렬로 표현하고, 역행렬을 이용해서 위와 같이 정리할 수 있습니다.

그리고, 역행렬을 구해버리면

끝나는 거지요. 여기서 몇몇 파라미터는 예전에 제가 사용한 것으로 고정하고

위 식에서

를 대입하여 정리하면

를 얻을 수 있습니다. 여기서 실제로는 저렇게 수치화시키지 않습니다. MATLAB/Simulink에서 변수로 사용할 수 있도록 하는데요. 우리는 그저 연습이니, 빨리(^^) 구현할 수 있도록 그냥 숫자로 모두 표현해버린 것입니다.

Simulink로 동역학 표현 !

사람들 마다 사용하는 목적이 다르겠지만, simulink는 이런 경우 미분방정식을 풀기 위해서 사용한다고 말할 수 있습니다. 그만큼 미분방정식의 형태로 표현되는 동역학 방정식의 구현에 큰 장점을 보이는데요.

이전의 시뮬링크기초에 관한 글(참조1, 참조2)을 참조하셔서 위의 형태로 꾸며주면 됩니다. 위에 표시된 부분은 제어입력 F인데, 지금은 제어입력을 구현하지 않으므로 '0'을 인가하면 됩니다.

위에 표시된 부분은 ddot_x인데요. 동역학과 잘 비교하시면 됩니다.

그래서 위와 같이 넣어주시고요.

위에 표시된 부분은 ddot_theta로서

로 구현하시면 됩니다.

저부분은 theta의 초기치를 넣는 부분입니다. 초기치를 주지 않으면, 시뮬레이션 하면 가만히 '0'의 상태에 머물러 있는것을 볼 수 있거든요. 그래서 좀 움직이라고, 초기치를 줘야죠.

simulink내부는 라디안단위이기 때문에 10도의 초기치를 위와같이 표현해서 넣어둡니다.

그리고, 위 gain은 라디안으로 각도를 관찰해야해서 좀 불편하니까 Scope에서는 도의 단위로 확인할 수 있도록 변환하는 것입니다.

그리고, 시뮬레이션을 시작해서 결과를 확인하면

이렇네요. 그냥 진자를 어떤 각도에 두고 자유운동을 시켜본겁니다.^^


댓글을 달아 주세요

  1. airman 2010.06.19 02:08 신고

    안녕하세요 요즘 자동제어 공부를 하고 있는데 정말로 많은 도움 받고 있습니다 ㅎㅎ
    그런데 한가지 궁금한게... 미방을 시뮬링크 말고 직접 풀순 없나요?
    룬지 쿠타 알고리즘이나 매틀랩 내장함수를 이용해서요

    • BlogIcon PinkWink 2010.06.19 13:52 신고

      네.. MATLAB 자체적으로 여러함수를 제공하고 있습니다. MATLAB document (command window 에서 doc이라고 입력하시면 됩니다.) 에서 ode45라고 검색해보시면 ode45뿐아니라 관련함수와 사용법이 나타납니다.^^

  2. AmOz 2010.12.30 15:08 신고

    위의 정리된 식들을 주욱 따라가다 보니, 역행렬 구하는 과정에서 계산 오류가 있는 것 같습니다. 1/(ad-bc) 에서 bc가 (m^2)(l^2) 인데, 따라서
    1/{(M+m)(ml^2+I)-(m^2)(l^2)} 가 되어야 할것입니다.

    • BlogIcon PinkWink 2010.12.31 16:46 신고

      네.. 말씀하신데로 부호에 오류가 보이네요.ㅠㅠ
      다행이 상수값부분인데다.. 원래 안정한 시스템이다 보니 제어기는 그냥 동작했던 모양이네요...
      하여간 발견해서 알려주셔서 감사합니다.^^

  3. 정인 2011.01.12 05:49 신고

    200x200x200 크기의 매트릭스를 z축에대해서 100도 회전하고 x축에관해서 10회전시키고 싶은 경우에 이용할 수 있는 매트랩 내장함수나 다른 방법이있는지 알고 싶어서요... 부탁^^;;

    • BlogIcon PinkWink 2011.01.12 13:52 신고

      일단... 200*200*200의 행렬의 회전이라는 것의 정의가 무엇인가요?ㅠㅠ
      보통 회전행렬은 X,Y,Z축 좌표(벡터로 표현된 그러니까 3*1짜리 행렬이죠)의 회전을 계산할때 사용하는 것으로 알고 있습니다만... 행렬이 회전한다는 개념을 모르겠습니다. ㅠㅠ

  4. 정인 2011.01.13 03:55 신고

    앗..질문이 구체적이 못해서 죄송합니다. Matrix라는 말을 행렬이라고 표현했는데요. 제코드가 일단 텍스트화일로부터 데이터를 읽어들이는데, 그 데이타크기가 (x=200,y=200,z=200)이고, 데이타들이 x,y,z각각의 위치에 저장되어있습니다. 이 데이터를 meshgrid를 사용하면 3D그래픽화가능한데, 원점이 (0,0,0)인데요, 저는 회전축을을 (x=100,y=100,z=100)으로 이동시킨상태에서 z축에대해 100도회전하고, x축에관해 10도 회전시켜보고 싶어서요... 이렇게하려면 어떤 방법이있을지 여쭤보고싶습니다.

    • BlogIcon PinkWink 2011.01.13 05:36 신고

      만약 해당 데이터를 X,Y,Z 순서로 1열 2열 3열에 배치해서 200*3크기의 행렬을 만들었다고 치고 그 이름을 Data라고 하죠. Z축중심의 회전행렬을 RotZ라 하고, X축 중심의 회전행렬을 RotX라고 하죠. (회전행렬의 구성이야 검색하시면 아실테니)

      Data*RotZ*RotX이라고 하면.. 간단히 구현될듯합니다만...
      (결과도 200*3의 크기를 가지구요)

      현재 새벽에 테스트없이 적은 댓글이니 바로될지는 잘 모르겠습니다만...

  5. 정인 2011.01.13 07:41 신고

    예. 그렇게 해보았는데, array index가 integer가 아니라는 오류때문에...

    • BlogIcon PinkWink 2011.01.13 10:30 신고

      200개까지는 못해보고.. 대략 10*3정도 크기의 Data행렬을 임의로 만들어서 테스트해보았는데, 문제가 없네요. 에러요인은 저에게 말씀하시지 않은 다른 곳에 있을 수도 있습니다. array index라는 것은 행렬내에서 위치를 의미합니다. 그래서 반듯이 정수여야하는 것인데요.

  6. 정인 2011.01.13 11:27 신고

    저는 잘 안되서 그러는데요. 계속 에러나오고요..
    ??? Subscript indices must either be real positive integers or logicals.

    Error in ==> test at 31
    val_voxel(updX, updY, updZ)=temp_val_dose(ix,jy,kz);
    라구요.
    제 코드좀 봐주시구요. 도움좀 주세요. 여기서 막혀서 진행이 안되서.. 도움을 바랍니다.

    count = 1;
    iteration = 200;
    for kz = 1:iteration
    for jy = 1:iteration
    for ix = 1:iteration
    val_test = count;
    temp_val_dose(ix,jy,kz) = val_test;
    count = count + 1;
    end
    end
    end

    the = 108*(2*pi)/360;
    phi = 10*(2*pi)/360;

    for kz = 1:iteration
    for jy = 1:iteration
    for ix = 1:iteration

    cnt = iteration;
    updX = ix*cos(the)*cos(phi)-jy*cos(the)*sin(phi)-kz*sin(the);
    updY = jy*cos(phi);
    updZ = ix*sin(the)*cos(phi)-jy*sin(the)*sin(phi)+kz*cos(the);

    val_voxel(updX, updY, updZ)=temp_val_dose(ix,jy,kz);
    cnt
    end
    end
    end

    • BlogIcon PinkWink 2011.01.13 16:42 신고

      제가 드릴 수 있는 도움은 코드전체를 봐드릴수있는 것이 아닙니다.
      그러나.. 살짝 코드를 봤을때...
      val_voxel(updX, updY, updZ)=temp_val_dose(ix,jy,kz);
      의 코드만 본다면...
      val_voxel(updX, updY, updZ)가 문제가 됩니다. 왜냐하면
      이 표현은 행렬에서 위치를 의미하는데
      updX, updY, updZ가 코사인과 사인의 계산으로 인해 정수가 아니기 때문입니다. 거기엔 몇번째라는 의미의 변수가 들어가야할것입니다.
      현재 서울로 이사준비(오늘밤에 떠나기때문에) MATLAB을 돌려볼수없습니다.

  7. 정인 2011.01.14 00:32 신고

    예.. 코드전체를 봐달란 뜻은 아니였는데, 핑크윙크님께선 된다고 하시니까 제건 왜 안되는지 알고 싶었습니다. 제 생각은 array index를 좌표상의 공간으로 보고 회전변환을 시키면 되지않을까했었거든요...

    • BlogIcon PinkWink 2011.01.15 15:18 신고

      이제 서울 이사를 마쳤네요...
      일단 아래에 starrynight님께서 좋은 소스를 공개해 두셨습니다.
      (일부 이모티콘으로 표현된것이 있지만 그 부분은 핵심이 되는 것은 아니기때문에 보시는데 괜찮으실겁니다.)
      starrynight님의 코드를 보시면 회전행렬(starrynight님의 회전행렬은 theta와 phi를 모두 회전하도록 두 회전행렬이 곱해져있습니다.)을 좌표와 곱하고 있습니다.
      제가 많이 도와드리지 못해 먼저 죄송합니다.ㅠㅠ

  8. BlogIcon starrynight 2011.01.14 16:34 신고

    정인 //

    예전에 만들었던 코드인데, 지금 하고자 하시는 것과 비슷한 것 같아서 올려보입니다. 혹시나 도움 되실까 하여..

    clear all; close all; clc;
    d = randn(3,200);

    theta_ = deg2rad(0);
    phi_ = deg2rad(1);
    psi_ = deg2rad(-2.5);

    rotMat = [...
    cos(theta_) * cos(psi_),...
    -cos(phi_)*sin(psi_)+sin(phi_)*sin(theta_)*cos(psi_),...
    sin(phi_)*sin(psi_)+cos(phi_)*sin(theta_)*cos(psi_);
    cos(theta_)*sin(psi_),...
    cos(phi_)*cos(psi_)+sin(phi_)*sin(theta_)*sin(psi_),...
    -sin(phi_)*cos(psi_)+cos(phi_)*sin(theta_)*sin(psi_);
    -sin(theta_),...
    sin(phi_)*cos(theta_),...
    cos(phi_)*cos(theta_);];

    D = rotMat*d;

    figure('color','w'),
    plot3(d(1,:),d(2,:),d(3,:),'r*'); box on; hold on; axis image; grid on;
    xlabel('X'); ylabel('Y'); zlabel('Z');
    camproj('perspective'); view(gca,[-50 30]);
    xlim([-max(d(:)),max(d(:))]);
    ylim([-max(d(:)),max(d(:))]);
    zlim([-max(d(:)),max(d(:))]);
    plot3(D(1,:),D(2,:),D(3,:),'b^')

    (핑크윙크님께서 기분 상해하시면 어쩌나 조마조마하며 ㅎㅎ;;)

  9. 정인 2011.01.16 10:14 신고

    보여주신 코드는 감사합니다. 그런데도 제가 아직 잘 이해를 못해서... 귀찮으실지모르지만, 제경우에 맞게 코드를 다시변형해 주시면 안될까요?
    저는 단순히 3차원회전행렬변환을 z축변환 x축변환을 곱해서 다음과같은 변환행렬을 사용했는데요.
    [ cos(the)*cos(phi)-cos(the)*sin(phi)-sin(the);
    cos(phi);
    sin(the)*cos(phi)-sin(the)*sin(phi)+cos(the); ];
    그리고, 200x200x200인 matrix를 z축으로 100도 x축으로 10도 회전시키는 경우일때를 가정해서, 회전시켜보고싶은데요. 저 코드를 어떻게 하면좋을지 도움을 부탁드립니다.

    • BlogIcon PinkWink 2011.01.16 14:24 신고

      p1 = [0 0 0];
      p2 = [1 0 0];
      p3 = [0 1 0];
      p4 = [0 0 1];
      p5 = [1 1 1];
      p6 = [1 1 0];

      points = [p1; p2; p3; p4; p5; p6];

      th_x = 10*pi/180;
      th_z = 45*pi/180;

      rotX = [1 0 0; 0 cos(th_x) -sin(th_x); 0 sin(th_x) cos(th_x)];
      rotZ = [cos(th_z) -sin(th_z) 0; sin(th_z) cos(th_z) 0; 0 0 1];

      NewPoints = points*rotX*rotZ

      그냥 간단히 해보았습니다. 좌표는 6개정도가 있다고 했고 (개수야 늘리시면 되니까요. 너무 많으면 for문을 사용하시면 될테니) 그걸 그냥 하나의 행렬로 합치고, 그러면 저는 6*3의 크기가 되고 말그대로 회전행렬과 곱했습니다. 그리고 마지막 결과가 newPoints라는 변수에 들어가있는데 혹시 그대로 사용하면 좋겠지만 빼내야한다면 역시 반복문을 이용해 첫줄은 newP1 이렇게 빼내면 될것입니다.

      그리고, 말씀하신 코드에는 적용할 수가 없었습니다. x축 z축의 회전행렬을 곱해도 3*3의 크기를 가지게 되는데, 지금 보여주신 코드는 3*1의 크기를 가지는 듯 해서 수정하기가 좀 어려웠습니다.

      그리고, 댓글 첫줄에 "귀찮으실지모르지만"이라고 하셨는데 안귀찮습니다. 다만 '정인'님께서 원하시는걸 정확히 알려드리지 못해 좀 안타까울 뿐입니다. 추운날씨 감기조심하세요^^

  10. 정인 2011.01.17 09:02 신고

    따뜻하신 말씀에 감사합니다. 보여주신 코드를 실험해보았는데요. 핑크윙크님과 별밤님께서는 점(1,0,0)을 X,Z축 변환해서 점(0.707,-0.696,0.122)으로 이동하게되는데요. 거기까지는 이해되는데요. 제가 원하는 것은 점(1,1,1)에 어떤 데이타 값이 들어있습니다. 4차원데이타라고 부르던데... 해서 점(1,1,1)은 변환해서 위치가점(0.7,-0.6,0.1)로 변환된 위치에 점(1,1,1)에 들어있는 데이타값을 변화시키지않고 위치 시키려고하는데, array index가 정수만 가능한터라 데이터값을 점(0.7,-0.6,0.1)에 매핑하고 싶은데, 불가능한질문인가요?

  11. 2011.01.17 09:18

    비밀댓글입니다

    • BlogIcon PinkWink 2011.01.17 09:32 신고

      간단합니다. 다만 지금은 업무중이라 오늘 밤 (혹시 오늘 밤에 회식이 있다면 내일 밤에 말씀드리겠습니다.)에 댓글로는 안되어서 포스팅을 할려고 합니다.^^

  12. 정인 2011.01.18 16:06 신고

    네.. 감사합니다.

  13. 불타는가슴 2011.05.16 12:50 신고

    실행하고 나서 scope를 더블클릭 하면 마지막 캡쳐 부분의 그림이 뜨잖아요
    근데 제가 그 그림의 axis를 조정하려고 그 그림에서 마우스 오른쪽클릭을
    하니까 Axes properties 라구 조절하는 부분이 있더라구요. 근데 여기선 Y축만 조절이 되고 X축은 없던데... X축까지 axis를 조절 하려면 어떻게 해야되나요??

    • BlogIcon PinkWink 2011.05.16 12:58 신고

      네.. 안됩니다.
      X축은 그저 시뮬레이션 시간으로만 조절이 가능합니다.

    • 불타는가슴 2011.05.16 16:40 신고

      아.. 그렇군요..
      제가 sumulink는 이번이 처음이라 모르는게 많습니다..

      일단 simulation 시간을 조절해 보려구 실행 아이콘

      바로옆에 simulation stop time을 조절해 보았는데

      변화가 없었는데... 그렇다면 simulation 시간을 어떻게 하나요???

    • BlogIcon PinkWink 2011.05.17 07:53 신고

      시뮬링크화면 상단 플레이버튼 옆에 숫자가 있는 곳을 수정하면됩니다.

  14. 요니아빠 2012.07.13 01:30 신고

    안녕하세요 정말 잘 배우고 있습니다.
    질문이 있어요, 힘은 '0'이고, 일정 각도(초기치)에 자유진동을 하게끔 설정되었는데
    결과 그래프를 보면, 왜 변위 'x'도 저렇게 진동하게 되는걸까요?
    진자의 진동이 미세하게나마 카트에 진동을 일으킨다고 해석해야되나요?

    • BlogIcon PinkWink 2012.07.17 08:40 신고

      실제와 흡사하게 시뮬레이션이 되었기 때문입니다.
      상상만 해봐도, 진자가 흔들거릴때, 그 진동이 카트(혹은 암에 전달될거니까요)

  15. kmh8075 2012.11.17 15:02 신고

    안녕하세요 학교에서 시뮬링크를 이용한 프로젝트를 수행중에 꼭 필요한 부분이 생겨 질문드립니다.

    혹시 input을 반복되는 step function으로 설정하는 방법 아시나요?

    0에서 1초까지 3
    1에서 2초까지 0
    2에서 3초까지 3
    3에서 4초까지 0
    이런 식으로 무한 반복되는 input을 입력하는 방법을 아신다며 좀 알려주시겠어요?

    • BlogIcon PinkWink 2012.11.18 16:02 신고

      시뮬링크에서 직접 혹은 mfile로해서 함수로 불러서든. 직접 반복되는 입력을 만드셔야할듯합니다. 혹은 제 블로그에서도 한번 소개한적이 있는 입력신호편집기를 사용해보시는것도 좋을듯하네요. 이름이 아마 signal builder였지 않나 하는데요

  16. sillf 2014.12.20 18:36 신고

    안녕하세요? 진자운동을 시뮬링크로 구성하는데요

    저위의 Fcn 과 먹스가 차원이 다르다고...output port 1 of mux is one dimensional vector with 3 elements 라고 뜨는데 뭐가 문제인지 잘모르겠습니다.

    • BlogIcon PinkWink 2014.12.21 21:53 신고

      mux에 연결된 선로가 가지는 데이터의 차원이 맞지 않는 모양이군요... 이것 말고는 드릴수 있는 말이 없습니다.ㅠㅠ

  17. 승연 2017.05.26 15:02 신고

    안녕하세요! 제어시스템공학 공부중인 학생입니다.

    위 시뮬링크는 선형 시스템의 시뮬레이션인데요

    똑같은 경우의 비선형 시스템 시뮬링크 구성은 어떻게 될까요??

    • BlogIcon PinkWink 2017.05.26 17:42 신고

      흠...ㅠㅠ.
      아마.. 비선형 모델이 있다면 (수학적으로표현된) 그 모델을 시뮬링크로 꾸미면 되겠지요ㅠㅠ.