루덴스코드 Blog

코딩과 교육/전문코딩 +51

소수구하는 알고리즘


예전에 파이썬으로 했던 소수 알고리즘이다.

2개의 수를 넣으면 그 두 수 사이의 모든 소수를 찾아서 화면에 출력한다.

실제로 소수를 찾는 시간은 매우 빠르지만 너무 많은 수를 출력하게 되면 화면 출력시간때문에 지연시간이 길어진다.



  • 2로 나누어 떨어지는 수는 소수가 아니므로 통과

  • 3으로 나누어 떨어지는 수는 소수가 아니므로 통과

  • 그러므로 3 이상의 모든 소수는 6k-1 또는 6k+1 에 해당한다. 모든 수에 대해서 % 연산을 수행하지 않고 6k-1, 6k+1 에 대해서만 수행한다.

  • 1부터 자기자신(N)까지의 모든 수를 대상으로 "N%(1...N)" 을 수행하지 않고, sqrt(N) 까지만 수행한다. 약수가 있다면 하나는 작은 수이고 다른 하나는 큰 수일텐데, 작은 수는 아무리 커도 sqrt(N) 보다는 작다. sqrt(N)*sqrt(N) = N 이므로 sqrt(N) 까지만 반복한다.

위의 조건에 따라 % 연산 결과가 0 이 되는 경우가 하나라도 발생하면 소수가 아니므로 그 상태에서 break; 를 수행


위의 조건을 따라 프로그램하면 다음과 같은 코드가 된다. (자바 코드)


package com.javalex.day3_1_primenumber;


public class PrimeNumberMake2 {

public static void primeMake(long startNumber, long endNumber) {

if (startNumber<=2 && endNumber >= 2) System.out.println(2 + "는 소수입니다");

if (startNumber<=3 && endNumber >= 3) System.out.println(3 + "는 소수입니다");

long modifiedNumber = 5;

if (startNumber > 5 ) modifiedNumber = startNumber;

for (long i = modifiedNumber; i<=endNumber; i++) {

boolean priNum = true;

if ((i%2==0)||(i%3==0)) priNum = false;

else {

for (int j=6; j<=(int)Math.sqrt(i)+1; j=j+6 ) {

if ((i%(j-1)==0)||(i%(j+1)==0)) { priNum = false; break; }

}

}

if (priNum==true) {

System.out.println(i+ "는 소수입니다");

}

}

}

}



아래는 위의 알고리즘이 적용되지 않은 소수 구하는 코드다. N 이 소수인것을 알기 위해 2 부터 N-1 까지 모두 나눠보고 0 으로 나누어서 떨어지는지를 검사한다. 보통 for 문을 배우면서 소수를 구하는 프로그램을 짜라고 시켜보면 이런식의 프로그램을 짠다. 


package com.javalex.day3_1_primenumber;


public class PrimeNumberMake0 {

public static void primeMake(long startNumber, long endNumber) {


if (startNumber<=2 && endNumber >= 2) System.out.println(2 + "는 소수입니다");

if (startNumber<=3 && endNumber >= 3) System.out.println(3 + "는 소수입니다");


long modifiedNumber = 4;

if (startNumber > 4 ) modifiedNumber = startNumber;

for (long i = modifiedNumber ; i<endNumber; i++) {

int priNum = 0;

for (int j=2; j<i; j++ ) {

if (i%j==0) priNum++;

}

if (priNum==0) System.out.println(i+ "는 소수입니다");

}

}

}



이 두개의 코드를 돌려서 걸리는 시간을 확인해 보았다. 이때 사용한 소스는 다음과 같다.


package com.javalex.day3_1_primenumber;


public class MainPrime {


public static void main(String[] args) {

long sNum = 100_000_000L;

long eNum = 100_000_100L;

long startTime;

startTime = System.currentTimeMillis();

PrimeNumberMake0.primeMake(sNum, eNum);

System.out.println((System.currentTimeMillis()-startTime)+"ms");


startTime = System.currentTimeMillis();

PrimeNumberMake2.primeMake(sNum, eNum);

System.out.println((System.currentTimeMillis()-startTime)+"ms");

}

}



결과는 다음과 같다.


100000007는 소수입니다

100000037는 소수입니다

100000039는 소수입니다

100000049는 소수입니다

100000073는 소수입니다

100000081는 소수입니다

194011ms

100000007는 소수입니다

100000037는 소수입니다

100000039는 소수입니다

100000049는 소수입니다

100000073는 소수입니다

100000081는 소수입니다

3ms


아무런 생각없이 그냥 짠 코드의 실행시간은 194초가 걸렸다. 내 PC 의 CPU 가 i3 라서 더 느리게 나왔다. 아마도 PC 성능이 좋으면 이보다는 조금 빨리 실행이 될 것이다. 반면 2의 배수, 3의 배수를 제외하고 6k+-1 의 항목만 검사하면서 2부터 N 까지가 아닌 까지 검사하면고, 소수가 아니면 바로 break 를 써서 빠진 결과는 3ms 가 나왔다.


194초 VS 0.003


저 차이가 바로 프로그램에 수학과 생각이 들어가야 하는 이유다.

Comment +0

아두이노를 이용한 PID제어기설계 : 

2013_Arduino PID Lab_0.pdf



아도이노 이용한 Ball 밸런스 : 

Balancing of a Ball on Beam using Arduino as a PID controller.hw



PID 제어 아두이노 사용 : 

PID-Control-Ardunio.pdf



PID제어 1, 2 : 

PID제어1.hwp

PID제어2.hwp





Comment +0

2014년 한양대에서 강의하는 확률과 통계 강의노트입니다.

 
필요한 부분을 출력해서 노트로 만들어 오시면 됩니다. 자료는 PDF 로 만들어져있으면 pdf view 프로그램을 통해 읽거나 프린트할 수 있습니다. 만약 pdf 를 읽는 방법을 모른다면 다음 링크로 가서 프로그램을 설치하시기 바랍니다.

http://www.foxitsoftware.com/kr/

프로그램 설치 후 위에서 다운받을 파일을 읽고, 프린트할 수 있습니다. 


강의시간에 늦지 않도록 노력해 주시기 바랍니다. 


 

Comment +0

2013년도 확률과 통계 강의 노트입니다.


기말고사 시간과 장소 공지

### 2013.06.19.수요일 오전 9시 - 11시 , 제1학술관 102호

### 계산기 지참할 것 - * 휴대폰의 계산기 프로그램 사용 불가, * 시험중인 타인의 계산기 빌릴 수 없음

### 06.24.월 오전 9-11시 사이 1공학관 5-315 에서 성적 상담, 상담전 반, 학번, 이름, 이상한 부분을 메일로 통해 알릴것 ( calltimy@hanyang.ac.kr )



 

 

확률과통계_01.pdf


 

 

확률과통계_02.pdf


 

 

확률과통계_03.pdf


 

 

확률과통계_04.pdf


 

 

확률과통계_05.pdf


 

 

확률과통계_06.pdf


 

 

확률과통계_07.pdf


 

 

확률과통계_08.pdf


 

 

확률과통계_09.pdf


 

 

확률과통계_10.pdf


 

 

확률과통계_11.pdf


 

 

확률과통계_12.pdf


 

 

확률과통계_13.pdf


 

 

확률과통계_14.pdf


 

 

확률과통계_15.pdf


 

 

확률과통계_16.pdf


 

 

확률과통계_17.pdf


 

 

확률과통계_18.pdf


 

 

확률과통계_19.pdf


  
  
  
  

한학기 동안 수고했습니다. 익숙치 않은 영어책으로 공부하는게 어려웠으리라 생각됩니다. 가능한 쉽게 풀어보려고 별도의 강의안을 만들어 봤습니다. 확통강의는 여러해 했지만 강의안을 프린트물로 정리한 것은 처음입니다. 도움이 되었기를 바랍니다.


무더운 여름 건강하게 지내시길 바랍니다.



Comment +0

확률과 통계 기말고사를 6월 15일 수요일 오후 7시부터 8시30분까지 시행하겠습니다.
장소는 아직 결정되지 않았으며 강의 시간에 알려드리도록 하겠습니다.


Comment +0

확률과 통계 시험 시간을 다음과 같이 알려드립니다. 장소는 추후 공지합니다.


23일 토요일 오전 9시-11시


iPhone 에서 작성된 글입니다.

Comment +0

라이브러리 파일입니다.


Comment +0

PDF Viewer Link : Foxit Reader 를 쓰면 편하다.
http://www.foxitsoftware.com/pdf/reader/

데이터북 전체를 한번에 보기


(1) 74LS 시리즈 (TTL)
(2) 74HC/HCT 시리즈 (CMOS/TTL)
(3) 74AC/ACT 시리즈 (CMOS/TTL)
(4) 4000/4500 시리즈 (CMOS)


Comment +0



강의평가시 사용할 자료입니다.
개별적인 설문도 조금 포함되어 있습니다.

Comment +0

실기평가용 PADS 라이브러리입니다.

Comment +0

PADS 로 만드는 8051 회로도 과제 입니다.



Comment +0

이 파일을 받으시고 참고하시기 바랍니다.

Comment +0

첨부된 파일을 참고 하세요.

Comment +1

호서대학교에서 공학입문설계에 사용하는 자료입니다.
PPT 로 만들어진 자료를 올립니다.
참고 하시기 바랍니다.
2009년도 1학기에 사용하는 내용입니다.
01 번은 강의 첫날 알려준 수업관련 성적, 배점, 출결 등에 관한 내용이었습니다. 여기에 올리지 않습니다.


PowerPoint 2007 이 설치되지 않은 경우 PowerPointer Viewer 를 설치하면 파일의 내용을 볼수 있습니다. Viewer 는 MS에서 무료로 제공해주고 있습니다. [LINK]

Comment +0

안산일대 2009년 CAD II 실기평가 자료
사용자 삽입 이미지

Comment +0

안산일대 중간고사 시험내용입니다. 2008년 2학기 자료입니다.



2008년도 2학기 C언어 중간고사 시험지입니다. 참고하여 시험준비하세요.

Comment +0

[START]

USB 를 이용한 기본 실험환경를 만들어 본다.

PC 에서 USB 를 사용해 외부 소형 보드에 연결하고, 보드에 붙어있는 LED 를 켜고 끄는 작업을 진행한다.


[참고사이트]

  1. FTDI 사이트
    1. http://www.ftdichip.com/
    2. USB 사용할 수 있는 칩제조사
    3. Evaluation Board 자료가 있다.
    4. FT245R, FT245BM 에 대한 자료
    5. http://www.ftdichip.com/Drivers/VCP.htm
    6. VCP (Virtual Com Port) Driver URL, D2XX Drivers URL,
  2. DLP Design
    1. http://www.dlpdesign.com/
    2. Visual C++ 사용한 소스코드제공
    3. http://www.dlpdesign.com/usb/usb245.shtml
  3. M2CV (국내회사)
    1. http://www.m2cv.co.kr/
    2. FT245 USB Module
    3. DLP Design 의 DLP-USB232M-G 과 동일
  4. 디바이스마트
    1. http://www.devicemart.co.kr/
    2. FT245 로 검색하면 관련 칩과 모듈이 검색됨

결론 : (잠정적)

회로도와 부품들을 살펴보고 PCB를 제작, 직접 모듈을 제작할 수도 있으나 부품을 구하기가 번거롭고, 제작시 칩형 부품들의 납땜에 에러가 발생할 것이 우려되어 모듈을 직접 구입하여 동작시켜 보기로 한다.


M2CV









DLP Design



FTDI














USB 를 이용한 제어 실험 - 첫번째
http://electoy.tistory.com/132
JelicleLim(2008.10.10)

Comment +0

사용자 삽입 이미지

엑손(exon)은 DNA 가닥상 단백질의 발현정보를 code 하고 있는 부분이고, 인트론(intron) 은 아직까지는 정확한 해석이 내려지지 않고 있는 부분이다. 일반적으로는 exon 부분을 보호하기 위한 일종의 껍질 같은 역할을 한다고 여겨진다.
인간 DNA 의 경우 염기서열중 5660만개 염기가 단백질을 합성하는 유전자 부분이고, 이는 전체 28억 3천만개의 유전자를 구성하는 염기서열 중 2%에 해당한다.(1)  
사용자 삽입 이미지

사용자 삽입 이미지

위 그림은 DNA 에 있는 정보를 mRNA가 가지고 오는 과정을 보여준다. 기본적으로 DNA의 정보가 mRNA로 전사되면서 모든 정보가 복사된다. 여기엔 불필요한 정보까지 복사되고 이 불필요한 정보를 스플라이싱과정을 거쳐 필요한 정보만으로 구성된 mRNA로 다시 만들게 된다.
사용자 삽입 이미지
이렇게 만들어진 mRNA 는 핵막을 통과할 수 있다.

사용자 삽입 이미지
이 mRNA 의 보다 정밀한 모습의 위와 같다. 기본적으로 엑손이 모인 중간 부분과 캡을 쓴 머리 부분, 그리고 여러개의 A 로 이어진 꼬리 부분이 있다. 여기서 모자같은 머리부분과 꼬리 부분을 제거하면 DNA의 정보를 그대로 담고 있는 부분이 된다.

Exon 과 Intron
http://electoy.tistory.com/127
JelicleLim(2008.9.4.)

P.S. (1) '유전자' -> '유전자를 구성하는 염기서열' 로 수정함, 댓글로 알려주셔서 고맙습니다. 

Comment +6

엔트로피(Entrophy)란 무엇인가?

사용자 삽입 이미지

에너지를 생각할 때 기본적인 법칙이 두 가지 있다. 그 첫째는 에너지란 만들어 질 수 없고, 없어지지 않는 다는 것이다. 열역학 제 1 법칙, 혹은 에너지 보존의 법칙이라고 불려지는 것인데, 모든 에너지는 항상 보존된다는 것을 의미한다.

에너지는 빛, 전기,그리고 열이라는 형태로 변하지만 그 양은 전체적으로 항상 일정하게 유지된다. 수력발전을 통해 물이 가지고 있던 위치에너지는 전기에너지로 변환되고, 그 전기에너지는 가정으로 전송되어 밤을 환하게 비추는 빛에너지가 되거나 방을 따뜻하게 만드는 난방장치를 통해 열어너지로 전환된다.

이렇게 전환된 에너지는 이전의 에너지와 이후의 에너지의 총량은 일정하다는 뜻이다. 물의 위치에너지는 100% 전기에너지로 전환되지 않는다. 일부는 전환되는 과정 중 마찰열을 발생시키고 물이 수증기가 되게 하는 에너지로도 사용된다. 즉 위치에너지는 전기에너지와 그 밖에 다른 에너지로 전환되고, 이 전환된 전기에너지는 물의 다시 원위치로 되돌리기에는 약간(혹은 그보다 상당히) 부족하게 된다는 것을 의미한다. 그리고 이후의 에너지는 이전과 같이 100% 사용 가능한 에너지가 아님을 의미한다. 에너지의 양은 같지만 사용 가능한 에너지의 양이 적어지는 것을 설명하는 것이 열역학 제 2 법칙, 또는 엔트로피 증가의 법칙이다.

[BIO] 엔트로피(Entrophy)란 무엇인가?
http://electoy.tistory.com/128
JelicleLim(2008.9.3)

Comment +0

2008년도 2학기 C언어 강의자료입니다.




교재 : "처음 만나는 C 프로그래밍 , Playing with C"
저자 : 우균, 창병모
출판사 : 교보문고

PLAYING WITH C 상세보기
우균 지음 | 교보문고 펴냄
C 입문서. 이 책은 C 언어의 입문과 상수, 입출력과 연산자, 제어구조와 함수, 배열, 포인터와 문자열 등의 내용을 다양한 예제를 통해 익힐 수 있도록 구성했다.

Comment +0

벡터의 삼중곱에 대한 내용이다. 내적과 외적을 사용한 삼중곱의 경우는 다음과 같다.

사용자 삽입 이미지

Comment +1

  • 궁금합니다. 2011.06.21 16:46 신고

    1.1 순환성이 왜 성립하는지 이해가 안되네요 어떻게 풀어야 같게 나오나요?
    ksng8130@naver.com 메일좀 주시면감사하겠습니다.

벡터의 내적 inner product 에 대한 기본 정의
사용자 삽입 이미지

벡터 내적의 물리적 의미는 상당히 심오하다. 고등학교 수학에서 간단하게 공식으로만 외웠던 위 한줄의 식이 의미하는 바를 제대로 이해하기 위해서는 상당히 내공이 필요하다. 내적이란 무엇인가? 어쩌면 철학적인 질문같기도 한 이 질문에 제대로 답할 수 있다면 그는 적어도 자신의 분야에서 제대로 된 내면의 진실을 조금은 접해본 사람이라고 할수도 있을 것이다.

벡터의 외적
사용자 삽입 이미지
외적은 내적과는 달리 크기와 방향을 함께 가지고 있는 벡터의 연산방법이다.

div A, 혹은 Del A, 라고 불리는 벡터의 발산은 다음과 같다.
사용자 삽입 이미지
벡터의 회전은 Del 혹은 Nabla 를 사용하여 나타낸다. Nabla 는 직교좌표 x, y, z에 관하여 각 성분이 ∂/∂x, ∂/∂y, ∂/∂z와 같은 형식적인 벡터로써 나타내어지는 미분연산자이다. W.R.해밀턴이 처음으로 사용한 연산자()이며, 기호는 ∇이다.
사용자 삽입 이미지

스톡스와 가우스는 각각 선적분을 면적분으로 면적분은 체적분으로 상호 변환하는 식을 만들었다. 이에 대해서는 별도의 공부가 필요할 듯...
사용자 삽입 이미지
사용자 삽입 이미지

벡터의 기본성질
http://electoy.tistory.com/123
JelicleLim(2008.8.16)

Comment +0

1961년 마셜 니렌버그와 조안 매테이는 어떤 종류의 RNA가 단백질을 만들도록 리보솜을 프로그램할 것이라는 가설을 세우고 실험을 진행했다. 결과적으로 그들은 우리딜산이라는 뉴클레오티드를 천연 RNA처럼 연결시킨 합성 RNA를 사용했을때 인공단백질을 생성하는 것을 발견하게 된다.
이들은 1965년까지 20가지 아미노산에 대해 각각의 세개 문자로 된 암호 61개를 모두 밝혀낸다. DNA의 4가지 뉴클레오티드는 ATGC 로 부르지만 RNA에서는 AUGC 로 부른다. DNA는 골격으로 디옥시리보스와 인산을 가지지만 RNA는 리보스와 인산의 골격을 갖는다.

사용자 삽입 이미지
표를 보면 알수 있듯이 CAU나 CAC나 같은 단백질을 만들어낸다. 즉, 4x4x4 라고 해서 64가지의 단백질이 만들어지는 것이 아니라 총 20가지의 단백질이 만들어지고, 세 종류의 코돈은 종결명령을 내리게 된다. AUG는 메티오닌을 만들거나 혹은 시작코돈, 개시코돈(Start Codon)의 역할을 한다.

생명과학이야기 상세보기
말론 &lt;b&gt;호클랜드&lt;/b&gt; 외 지음 | 진솔서적 펴냄
The way life works의 후속판에 해당하는 책으로 생명이 가지고 있는 본질적인 요소를 분명하고 알기쉽게 기술하였다. 생물학적인 예와 과학적인 방법을 많이 첨가하여 이해를 돕고 있다.


유전암호의 해독
http://electoy.tistory.com/122
JelicleLim(2008.8.13)

Comment +0

하나의 알파벳들이 모여 단어를 이루고, 그 단어가 모여 문장을 이루고, 그 문장이 모여 책을 이루고 책이 모여 전집을 이루듯이 생명체의 기본단위인 세포(유전체)는 가장 작은 단위의 뉴클레오티드가 모여 삼중체를 형성하고, 그 삼중체가 유전자를 형성하고, 그것이 다시 염색체를 형성하고 그것이 하나의 세포의 핵을 이루고 세포를 형성, 생명체를 구성하게 된다.

사용자 삽입 이미지
사용자 삽입 이미지
사용자 삽입 이미지
사용자 삽입 이미지
사용자 삽입 이미지

위와 같은 순서가 된다. 굳이 순서를 외우려고 할 필요는 없다. 다만 여기서 신경써야 할 것은 AT와 CG 의 4개의 뉴클레오티드의 결합이 정보를 만들게 되고 그 정보가 생명체의 가장 중요한 암호가 된다.
앞으로의 바이오인포매틱스에서 이 암호화된 정보를 어떻게 해석해 내느냐가 중요한 흐름이기도 하다.

황박사의 줄기세포같은 연구도 필요하겠지만 실상 DNA의 정보를 해석해내지 못하는 한 세포의 연구는 기껏해야 복제이상의 성과를 내지 못한다. 즉, 유전병이나 암의 치료제등은 아무리 줄기세포를 많이 가지고 있어도 별 소용이 없다는 말이다. 세포속 암호화된 DNA의 정보를 깨닫지 못한다면 말이다.

지금 하고 있는 연구는 이 DNA 정보를 분석하는 작업을 하고 있다. 지금은 초기단계라 간단한 실험정도만 하고 코돈(삼중체)의 해석을 하고 있지만 그것도 쉽지 않은 것이고, 딱히 주어진 텍스트가 없이 추정으로 진행해야 하는 부분이 많은 것도 사실이다.

뉴클레오티드에서 게놈까지
http://electoy.tistory.com/121
JelicleLim(2008.8.13)

Comment +0

17세기 식물은 성장은 흙에서 비롯된다고 사람들은 생각했다. 그리고 이에 대해 새로운 실험을 당시 플랑드르의 내과의사 잔 뱁티스타 반 헬몬트가 한다.
사용자 삽입 이미지

헬 몬트는 2kg 정도의 버드나무 가지를 90kg 의 흙에 심었다. 그리고 규칙적으로 물을 주면서 5년을 키운다. 나무가지의 무게는 75kg 이 된 반면 흙은 57g 만 줄어들었다. 여기는 헬몬트는 나무를 구성하는 물질이 흙에서 온 것이 아니라 물에서 온 것이라는 결론을 내린다.

헬몬트의 실험은 당시 사람들이 가지고 있던 잘못된 생각을 깨우치는 데 도움이 되었지만 그의 실험은 완벽했던 것은 아니다. 지금은 나무를 구성하는데 소요되는 것이 물과 흙 뿐만이 아니라 공기, 즉 탄소라는 것을 알게 되었지만 당시로는 그것까지는 미처 알수 없었던 것이다.


Comment +0

Ball and Beam problem using PID Control

In this digital control version of the ball and beam experiment, we are going to use the PID control method to design the digital controller. If you refer to the Ball and Beam Modeling page, the open-loop transfer function was derived as
사용자 삽입 이미지
.....m........mass of the ball...............0.11 kg
.....g........gravitational acceleration.....9.8 m/s^2
.....d........lever arm offset...............0.03 m
.....L........length of the beam.............1.0 m
.....R........radius of the ball.............0.015 m
.....J........ball's moment of inertia.......9.99e-6 kgm^2
.....R(s).....ball position coordinate (m)
.....theta(s).servo gear angle...............0.25 rad

The design criteria for this problem are:
  • Settling time less than 3 seconds
  • Overshoot less than 5%
Digital PID controller
If you refer to any of the PID control problem for continuous systems, the PID transfer function was expressed as
사용자 삽입 이미지
As you noticed the above transfer function was written in terms of s. For the digital PID control, we use the following transfer function in terms of z.
사용자 삽입 이미지
Discrete transfer function
The first thing to do here is to convert the above continuous system transfer function to discrete transfer function. To do this, we are going to use the Matlab function called c2dm. To use this c2dm, we need to specify four arguments: numerator and denominator matrices, sampling time (Ts), and the 'method'. You should already be familiar with how to enter numerator and denominator matrices. The sampling time should be smaller than 1/(30*BW) sec, where BW is the closed-loop bandwidth frequency. The method we will use is the zero-order hold ('zoh'). Assuming that the closed-loop bandwidth frequency is around 1 rad/sec, let the sampling time be 1/50 sec/sample. Now we are ready to use c2dm. Enter the following commands to an m-file.

.....m = 0.111;
.....R = 0.015;
.....g = -9.8;
.....L = 1.0;
.....d = 0.03;
.....J = 9.99e-6;

.....K = (m*g*d)/(L*(J/R^2+m));   %simplifies input

.....num = [-K];
.....den = [1 0 0];

.....Ts = 1/50;

.....[numDz,denDz]= c2dm (num,den,Ts,'zoh')

Running this m-file in the Matlab command window gives you the following matrices.

.....numDz =
..........1.0e-0.4 *
..........0    0.4200    0.4200

.....denDz =
..........1    -2    1

From these matrices, the discrete transfer function can be written as
사용자 삽입 이미지
Open-loop response
Now we will observe the ball's response to a step input of 0.25 m. To do this, enter the following commands to an new m-file and run it in the command window. You should see the following response.

.....numDz = 0.0001*[0.42 0.42];
.....denDz = [1 -2 1];

.....[x] = dstep(0.25*numDz,denDz,251);
.....t=0:0.02:5;
.....stairs(t,x)
사용자 삽입 이미지
From this plot, it is clear that the open-loop system is unstable causing the ball to roll off from the end of the beam.

Proportianal Control
Now we will add the proportional control (Kp) to the system and obtain the closed-loop system response. For now let Kp equal to 100 and see what happens to the response. Enter the following commands to an new m-file and run it in the command window.

.....numDz = 0.0001*[0.42 0.42];
.....denDz = [1 -2 1];

.....Kp=100;
.....[numDzC,denDzC]=cloop (Kp*numDz,denDz);

.....[x] = dstep(0.25*numDzC,denDzC,251);
.....t=0:0.02:5;
.....stairs(t,x)
사용자 삽입 이미지
As you can see, the addition of proportional control does not make the system stable. You may try to increase the proportional gain (Kp) and confirm that the system remains unstable.

Proportional-Derivative control
Now we will add a derivative term to the controller. Keep the proportional gain (Kp) equal to 100, and let the derivative gain (Kd) equal to 10. Copy the following code to an new m-file and run it to view the system response.

.....numDz = 0.0001*[0.42 0.42];
.....denDz = [1 -2 1];
.....Kp=100;
.....Kd=10;
.....numpd = [Kp+Kd -(Kp+2*Kd) Kd];
.....denpd = [1 1 0];
.....numDnew = conv(numDz,numpd);
.....denDnew = conv(denDz,denpd);
.....[numDnewC,denDnewC] = cloop(numDnew,denDnew);
.....[x] = dstep(0.25*numDnewC,denDnewC,251);
.....t=0:0.02:5;
.....stairs(t,x)
사용자 삽입 이미지
Now the system is stable, but the rise time is too long. From the PID Tutorial page, we see that the increasing the proportional gain (Kp) will decrease the rise time. Let's increase the proportional gain (Kp) to 1000 and see what happens. Change Kp in the above m-file from 100 to 1000 and rerun it in the command window. You should see the following step response.
사용자 삽입 이미지
As you can see, all of the design requirements are satisfied. For this particular problem, no implementation of an integral control was needed. But remember there is more than one solution for a control problem. For practice, you may try different P, I and D combinations to obtain a satisfactory response.

Comment +0

Ball & Beam Problem Using the State-space Design Method

The state-space representation of the ball and beam example is given below:
사용자 삽입 이미지
Remember, unlike the previous examples where we controlled the gear's angle to control the beam and ball, here we are controlling alpha-doubledot. By doing this we are essentially controlling a torque applied at the center of the beam by a motor. Therefore, we do not need a gear and lever system.

The design criteria for this problem are:
  • Settling time less than 3 seconds
  • Overshoot less than 5%
To see the derivation of the state-space equations for this problem refer to the ball and beam modeling page.

If you are interested in running an animation of this example based on the control techniques used in the state-space tutorial please go to the Ball & Beam Animation Page after completing this tutorial.

Full-State Feedback Controller

We will design a controller for this physical system that utilizes full-state feedback control. A schematic of this type of system is shown below:
사용자 삽입 이미지
Recall, that the characteristic polynomial for this closed-loop system is the determinant of (sI-(A-BK)), where s is the Laplace variable. For our system the A and B*K matrices are both 4x4. Hence, there should be four poles for our system. In designing our full-state feedback controller we can move these poles anywhere we want.

For our design we desire an overshoot of less than 5% which corresponds to a zeta of 0.7 (please refer to your textbook for the relationship between overshoot and damping ratio). On a root locus this criterion is represented as a 45 degree line emanating from the origin and extending out into the left-half plane. We want to place our poles on or beneath this line. Our next criterion is a settling time less than 3 seconds, which corresponds to a sigma = 4.6/Ts = 4.6/3 = 1.53, represented by a vertical line at -1.53 on the root locus. Anything beyond this line in the left-half plane is a suitable place for our poles. Therefore we will place our poles at -2+2i and -2-2i. We will place the other poles far to the left for now, so that they will not affect the response too much. To start with place them at -20 and -80. Now that we have our poles we can use Matlab to find the controller (K matrix) by using the place command. Copy the following code to an m-file to model the system and find the K matrix:

NOTE: Matlab commands from the control system toolbox are highlighted in red.


.....m = 0.111;
.....R = 0.015;
.....g = -9.8;
.....J = 9.99e-6;

.....H = -m*g/(J/(R^2)+m);
         
.....A=[0 1 0 0
........0 0 H 0
........0 0 0 1
........0 0 0 0];
.....B=[0;0;0;1];
.....C=[1 0 0 0];
.....D=[0];

.....p1=-2+2i;
.....p2=-2-2i;
.....p3=-20;
.....p4=-80;

.....K=place(A,B,[p1,p2,p3,p4])

Run your m-file and you should get the following output for the K matrix:

.....place: ndigits= 15

.....K =
..........1.0e+03 *
..........1.8286    1.0286    2.0080    0.1040

After adding the K matrix, the state space equations now become:
사용자 삽입 이미지
We can now simulate the closed-loop response to a 0.25m step input by using the lsim command. Add the following to your m-file:

.....T = 0:0.01:5;                  
.....U = 0.25*ones(size(T));             
.....[Y,X]=lsim(A-B*K,B,C,D,U,T);     
.....plot(T,Y)

Run your m-file and you should get the following plot:
사용자 삽입 이미지
From this plot we see that there is a large steady state error for which we will need to add reference input (explained in next section). However, the overshoot and settling time criteria are met. If we wanted to reduce the overshoot further than we would make the imaginary part of the pole smaller than the real part. Also, if we wanted a faster settling time we would move the poles further in the left-half plane. Feel free to experiment with the pole positions to see these trends.

Reference Input
Now we want to get rid of the steady-state error. In contrast to the other design methods, where we feedback the output and compare it to the reference input to compute an error, with a full-state feedback controller we are feeding back both states. We need to compute what the steady-state value of the states should be, multiply that by the chosen gain K, and use a new value as our reference for computing the input. This can be done by adding a constant gain Nbar after the reference. The schematic below shows this relationship:
사용자 삽입 이미지
Nbar can be found using the user-defined function rscale (copy it to the directory that your m-file is in). Copy the following to your m-file and run it to view the step response with Nbar added.

.....Nbar=rscale(A,B,C,D,K)

.....T = 0:0.01:5;                  
.....U = 0.25*ones(size(T));             
.....[Y,X]=lsim(A-B*K,B*Nbar,C,D,U,T);     
.....plot(T,Y)

Note: Non-standard Matlab commands used in this example are highlighted in green.

Your output should be:

.....place: ndigits= 15
.....Nbar =
..........1.8286e+03
사용자 삽입 이미지
Now the steady-state error is gone and all the design criteria are satisfied.

Note: A design problem does not necessarily have a unique answer. Using this method (or any other) may result in many different compensators. For practice you may want to go back and try to change the pole positions to see how the system responds.

Comment +0

Ball & Beam Problem Using Frequency Response Method

The open-loop transfer function of the plant for the ball and beam experiment is given below:
사용자 삽입 이미지
The design criteria for this problem are:
  • Settling time less than 3 seconds
  • Overshoot less than 5%
To see the derivation of the equations for this problem refer to the ball and beam modeling page. A schematic of the closed loop system with a controller is given below:
사용자 삽입 이미지

Open-loop Bode Plot
The main idea of frequency based design is to use the Bode plot of the open-loop transfer function to estimate the closed-loop response. Adding a controller to the system changes the open-loop Bode plot, therefore changing the closed-loop response. Let's first draw the bode plot for the original open-loop transfer function. Create an m-file with the following code and then run it in the Matlab command window:

.....m = 0.111;
.....R = 0.015;
.....g = -9.8;
.....L = 1.0;
.....d = 0.03;
.....J = 9.99e-6;

.....K = (m*g*d)/(L*(J/R^2+m));   %simplifies input

.....num = [-K];
.....den = [1 0 0];

.....bode(num,den)

NOTE: Matlab commands from the control system toolbox are highlighted in red.


You should get the following Bode plot:
사용자 삽입 이미지
From this plot we see that the phase margin is zero. Since the phase margin is defined as the change in open-loop phase shift necessary to make a closed-loop system stable this means that our zero phase margin indicates our system is unstable. We want to increase the phase margin and we can use a lead compensator controller to do this. For more information on Phase and Gain margins please refer to the Frequency Response Tutorial.

Phase-Lead Controller
A first order phase-lead compensator has the form given below:
사용자 삽입 이미지

The phase-lead compensator will add positive phase to our system over the frequency range 1/aT and 1/T, which are called the corner frequencies. The maximum added phase for one lead compensator is 90 degrees. For our controller design we need a percent overshoot of 5, which corresponds to a zeta of 0.7. Generally zeta * 100 will give you the minimum phase margin needed to obtain your desired overshoot. Therefore we require a phase margin greater than 70 degrees.

To obtain "T" and "a", the following steps can be used.

1. Determine the positive phase needed:

      We need at least 70 degrees from our controller.

2. Determine the frequency where the phase should be added (center frequency):

      In our case this is difficult to determine because the phase vs. frequency graph in the bode plot is a flat line. However, we have a relation between bandwidth frequency (wbw) and settling time (refer to the Bandwidth Frequency page for this equation) which tells us that wbw is approximately 1.92 rad/s. Therefore we want a center frequency just before this. For now we will choose 1.

3. Determine the constant "a" from the equation below, this determines the required space between the zero and the pole for the maximum phase added.
사용자 삽입 이미지
.....where phi refers to the desired phase margin. For 70 degrees, a = 0.0311.

4. Determine "T" and "aT" from the following equations:
사용자 삽입 이미지
사용자 삽입 이미지
..... For 70 degrees and center frequency (w) = 1, aT = 0.176 and T = 5.67

Now, we can add our lead controller to the system and view the bode plot. Remove the bode command from your m-file and add the following:

.....k=1;
.....numlead = k*[5.67 1];
.....denlead = [0.176 1];

.....numl = conv(num,numlead);
.....denl = conv(den,denlead);

.....bode(numl,denl)

You should get the following bode plot:
사용자 삽입 이미지

You can see that our phase margin is now 70 degrees. Let's check the closed-loop response to a step input of 0.25m. Add the following to your m-file:

.....[numcl,dencl] = cloop(numl,denl);
.....t=0:0.01:5;
.....step(0.25*numcl,dencl,t)

You should get the following plot:
사용자 삽입 이미지
Although the system is now stable and the overshoot is only slightly over 5%, the settling time is not satisfactory. Increasing the gain will increase the crossover frequency and make the response faster. Make k = 5, your response should look like:
사용자 삽입 이미지
사용자 삽입 이미지
The response is faster, however, the overshoot is much too high. Increasing the gain further will just make the overshoot worse.

Adding More Phase
We can increase our phase-lead compensator to decrease the overshoot. In order to make the iterative process easier use the following program. Create an m-file and copy the function from your web-browser into it (make sure the function command starts in the first column of the m-file).

.....function[ ] = phaseball()

.....%define TF
.....m = 0.111;
.....R = 0.015;
.....g = -9.8;
.....L = 1.0;
.....d = 0.03;
.....J = 9.99e-6;

.....K = (m*g*d)/(L*(J/R^2+m));   %simplifies input

.....num = [-K];
.....den = [1 0 0];

.....%ask user for controller information
.....pm = input('Phase Margin?.......');
.....w  = input('Center Frequency?...');
.....k  = input('Gain?...............');

.....%view compensated system bode plot
.....pmr = pm*pi/180;
.....a = (1 - sin(pmr))/(1+sin(pmr));
.....T = sqrt(a)/w;
.....aT = 1/(w*sqrt(a));

.....numlead = k*[aT 1];
.....denlead = [T 1];

.....numl=conv(num,numlead);
.....denl=conv(den,denlead);
.....figure
.....bode(numl,denl)

.....%view step response
.....[numcl,dencl]=cloop(numl,denl);
.....t=0:0.01:5;
.....figure
.....step(0.25*numcl,dencl,t)

With this m-file you can choose the phase margin, center frequency, and gain. Run your m-file with the following values and you should see the plots below on your screen.

.....Phase Margin?.......80
.....Center Frequency?...1
.....Gain?...............1
사용자 삽입 이미지
사용자 삽입 이미지
The overshoot is fine but the settling time is just a bit long. Try different numbers and see what happens.

Using the following values the design criteria was met.

.....Phase Margin?.......85
.....Center Frequency?...1.9
.....Gain?...............2
사용자 삽입 이미지
사용자 삽입 이미지
Note: A design problem does not necessarily have a unique answer. Using this method (or any other) may result in many different compensators. For practice you may want to go back and change the added phase, gain, or center frequency.

Comment +0

Solution to the Ball & Beam Problem Using Root Locus Method

The open-loop transfer function of the plant for the ball and beam experiment is given below:

사용자 삽입 이미지
The design criteria for this problem are:
  • Settling time less than 3 seconds
  • Overshoot less than 5%

To see the derivation of the equations for this problem refer to the ball and beam modeling page. A schematic of the closed loop system with a controller is given below:

사용자 삽입 이미지

Open-loop Root Locus

The main idea of the root locus design is to estimate the closed-loop response from the open-loop root locus plot. By adding zeroes and/or poles to the original system (adding a compensator), the root locus and thus the closed-loop response will be modified. Let us first view the root locus for the plant in open loop. Create an m-file with the following Matlab code in order to model the plant and plot the root locus.

.....m = 0.111;
.....R = 0.015;
.....g = -9.8;
.....L = 1.0;
.....d = 0.03;
.....J = 9.99e-6;

.....K = (m*g*d)/(L*(J/R^2+m));   %simplifies input

.....num = [-K];
.....den = [1 0 0];

.....rlocus(num,den)

NOTE: Matlab commands from the control system toolbox are highlighted in red.

Now, run the m-file and you should see the following root locus plot:

사용자 삽입 이미지
As you can see the system has two poles at the origin which go off to infinity along the imaginary axes.

The design criteria can also be plotted onto the root locus using the sgrid command. This command generates a grid of constant damping ratio and natural frequency. The damping ratio and natural frequency were found using the following equation, which relates the them to our percent overshoot (PO) and settling time (Ts) requirements:

사용자 삽입 이미지
사용자 삽입 이미지
Note, that the equation with Ts is found by assuming settled is when the response remains within 2% of it's final value. From these equation damping ratio and natural frequency were found to be 0.7 and 1.9 respectively.

.....sgrid(0.70, 1.9)
.....axis([-5 5 -2 2])

사용자 삽입 이미지
The area between the two dotted diagnol lines represents locations where the percent overshoot is less than 5%. The area outside the curved line represents locations where the setttling time is less than 3 seconds. Note that no region of the plot falls within the design criteria shown be these lines. To remedy this and bring the root locus into the left-hand plane for stability we will try adding a lead-compensator to the system.

Lead Controller
A first order lead compensator tends to shift the root locus into the left-hand plane. For a more detailed description of lead compensators refer to the Lead & Lag Compensator Design page. A lead compensator has the form given below:

사용자 삽입 이미지
where, the magnitude of zo is less than the magnitude of po.

Now, let us add the controller to the plant and view the root locus. We will position the zero near the origin to cancel out one of the poles. The pole of our compensator will be placed to the left of the origin to pull the root locus further into the left-hand plane. Add the following lines of Matlab code to your m-file.

.....zo = 0.01;
.....po = 5;

.....numlead = [1 zo];
.....denlead = [1 po];

.....numl = conv(num,numlead);
.....denl = conv(den,denlead);

.....rlocus(numl,denl)
.....sgrid(0.70, 1.9)


Run your m-file in the Matlab command window and you should see the following:

사용자 삽입 이미지
Now, the branches of the root locus are within our design criteria.

Selecting a Gain
Now that we have moved the root locus into the left-hand plane, we may select a gain that will satisfy our design requirements. We can use the rlocfind command to help us do this. Add the following onto the end of your m-file.

.....[kc,poles]=rlocfind(numl,denl)

Go to the plot and select a point near those indicated by the cross mark on the plot below:

사용자 삽입 이미지
You should see in the Matlab command window something similar to the following (your numbers will be slightly different):


.....selected_point =
.....  -2.4988+ 1.2493i
 
.....kc =
.....   37.3131
 
.....poles =
.....  -2.4950+ 1.2493i
.....  -2.4950- 1.2493i
.....  -0.0101       


Now, we can plot the response with this gain.

Plotting the Closed-loop Response
This value of kc can be put into the system and the closed-loop response to a step input of 0.25 m can be obtained. Add the following lines to your m-file to perform this analysis.

.....numl2 = kc*numl;
.....[numcl,dencl] = cloop(numl2,denl);
.....t=0:0.01:5;
.....figure
.....step(0.25*numcl,dencl,t)

Run your m-file and you select a point on the root locus similar to the selected point above. The step response should look like the following:

사용자 삽입 이미지
From this plot we see that when a 0.25m step input is given to the system both the settling time and percent overshoot design criteria are met.
Note: A design problem does not necessarily have a unique answer. Using this method (or any other) may result in many different compensators. Try running your m-file several more times selecting a different point each time and study the effect this has on the step response. For practice you may also want to go back to the original open-loop root locus and try to find other ways to add zeros and poles to get a better response.

Comment +0

Solution to the Ball & Beam Problem Using PID Control

The open-loop transfer function of the plant for the ball and beam experiment is given below:

사용자 삽입 이미지
The design criteria for this problem are:
  • Settling time less than 3 seconds
  • Overshoot less than 5%

To see the derivation of the equations for this problem refer to the ball and beam modeling page.


Closed-loop Representation
The block diagram for this example with a controller and unity feedback of the ball's position is shown below:

사용자 삽입 이미지
First, we will study the response of the system shown above when a proportional controller is used. Then, derivative and/or integral control will be added if necessary.

Recall, that the transfer function for a PID controller is:

사용자 삽입 이미지

Proportional Control
The closed-loop transfer function for proportional control with a proportional gain (kp) equal to 100, can be modeled by copying the following lines of Matlab code into an m-file (or a '.m' file located in the same directory as Matlab)

.....m = 0.111;
.....R = 0.015;
.....g = -9.8;
.....L = 1.0;
.....d = 0.03;
.....J = 9.99e-6;

.....K = (m*g*d)/(L*(J/R^2+m));   %simplifies input

.....num = [-K];
.....den = [1 0 0];

.....kp = 1;
.....numP = kp*num;

.....[numc, denc] = cloop(numP, den)

NOTE: Matlab commands from the control system toolbox are highlighted in red.

You numerator and denominator should be:

.....numc =..... 
.....         0         0    0.2100

.....denc =
.....    1.0000         0    0.2100

Now, we can model the system's response to a step input of 0.25 m. Add the following line of code to your m-file and run it:

.....step(0.25*numc,denc)

You should get the following output:

사용자 삽입 이미지
As, you can see the addition of proportional gain does not make the system stable. Try changing the value of kp and note that the system remains unstable.

Proportional-Derivative Control
Now, we will add a derivative term to the controller. Copy the following lines of code to an m-file and run it to view the system's response to this control method.

.....m = 0.111;
.....R = 0.015;
.....g = -9.8;
.....L = 1.0;
.....d = 0.03;
.....J = 9.99e-6;

.....K = (m*g*d)/(L*(J/R^2+m));   %simplifies input

.....num = [-K];
.....den = [1 0 0];

.....kp = 10;
.....kd = 10;
.....numPD = [kd kp];

.....numh = conv(num, numPD);
.....[numc, denc] = cloop(numh, den);
.....t=0:0.01:5;
.....step(0.25*numc,denc,t)

Your plot should be similar to the following:

사용자 삽입 이미지
Now the system is stable but the overshoot is much too high and the settling time needs to go down a bit. From the PID tutorial page in the section on characteristics of P, I, and D controllers, we see that by increasing kd we can lower overshoot and decrease the settling time slightly. Therefore, make kd = 20 in your m-file and run it again. Your output should be:
사용자 삽입 이미지
The overshoot criterion is met but the settling time needs to come down a bit. To decrease the settling time we may try increasing the kp slightly to increase the rise time. The derivative gain (kd) can also be increased to take off some of the overshoot that increasing kp will cause. After playing with the gains a bit, the following step response plot can be achieved with kp = 15 and kd = 40:
사용자 삽입 이미지
As you can see from the above plot all the control objectives have been met without the use of an integral controller (settling time for this example is considered achieved when the response is less than 2% of it's final value). Remember, that for a control problem there is more than one solution for the problem.

Comment +0