Otomatik Kontrol Ulusal Toplantısı, TOK2013, 26-28 Eylül 2013, Malatya AKILLI VE UYARLAMALI KONTROL SİSTEMLERİ 96 Otomatik Kontrol Ulusal Toplantısı, TOK2013, 26-28 Eylül 2013, Malatya Uyarlanabilir Sinirsel Bulanık Çıkarım Sistemi (ANFIS) ile Küçük Bir Turboprop Motor Üzerinde Yakıt Akışı Kontrol Sistemi Modellenmesi Işıl Yazar1, Emre Kıyak2, Fikret Çalışkan3, Kerim Kahraman4 1Mekatronik Bölümü Eskişehir Osmangazi Üniversitesi, Eskişehir [email protected] 2Havacılık Elektrik ve Elektroniği Bölümü Anadolu Üniversitesi, Eskişehir [email protected] 3Kontrol Mühendisliği Bölümü İstanbul Teknik Üniversitesi, İstanbul [email protected] 4TEI Tusaş Motor Sanayii A.Ş., Eskişehir [email protected] kompresör elemanının dönme hareketini sağlar. Türbinde Özetçe atmosfer basınç değerine kadar genleşen gaz kütlesi egzozdan Ülkemizde son dönemlerde savunma sanayii alanında ciddi dışarı atıldığında hava aracı için gerekli olan itkinin oluşması gelişmeler kaydedilmektedir. Gelişen teknolojiye bağlı olarak sağlanır. tasarlanmakta olan hava araçları da buna paralel olarak Turboprop motorlar literatürde iki başlık altında çeşitlilik göstermektedir. Özellikle İHA’lar (İnsansız Hava incelenmektedirler [2, 3]. Bunlar sırasıyla tek şaftlı turboprop Araçları) üzerinde ciddi çalışmalar yürütülmektedir. motorlar ve serbest türbinli turboprop motorlardır. Tek şaftlı Savunma, eğitim, gözlem gibi alanlarda kullanılan bu tip hava turboprop motorlarda pervane, kompresör-türbin yapısını da araçlarının itki ihtiyacını karşılamak için küçük tipte çeviren ortak bir şaft tarafından döndürülmektedir. Serbest turboprop motorlara gereksinim duyulmaktadır. Turboprop türbinli turboprop motorlarda ise halihazırdaki türbin motorlar, düşük ve orta irtifa için kullanımı ekonomik, sessiz yapısından ayrı olarak güç türbini olarak isimlendirilen bir çalışan gaz türbinli motorlardır. Çalışmada, TEI tarafından yapı pervane hareketinden sorumludur. Her iki yapıda da geliştirilen bir turboprop motorun verileri kullanılarak ANFIS türbin pervaneden daha hızlı döndüğü için pervanenin devir (Adaptive Neuro-Fuzzy Inference System) arayüzü sayısı şaft üzerinde bulunan dişli kutusu yardımı ile yardımıyla yakıt akışı kontrolü sağlanmıştır. Modelde hem ayarlanmaktadır. Şekil.1’de AIRBUS firmasına ait askeri geçici durum hem de kararlı durum süreçleri için yakıt akışı nakliye uçağı A400M’in TP400 Turboprop motoru modeli durumu incelenmiştir. görülmektedir [4]. 1. Giriş Günümüzde gaz türbinli motorların hem endüstri alanında hem de havacılık sektöründe geniş bir kullanım alanı bulunmaktadır. Endüstri alanındaki kullanımı ile enerji eldesi sağlanırken, havacılık sektöründeki kullanımı ile hava aracı için gereken itki sağlanmaktadır. Hava araçlarında kullanılan gaz türbinli motorlar şu şekilde sınıflandırılabilir: Turbojet, Turbofan, Ramjet, Pulsejet, Turboprop, Propfan [1]. Gaz türbinli motorlar kendi aralarında yapısal olarak farklılık gösterse de temel çalışma prensibi hepsinde aynıdır. Hava alığından içeri alınan havanın hızı düşürülerek basıncı bir miktar artırılır. Daha sonra kompresör içerisine geçen havanın sabit sıcaklıkta sıkıştırılma yoluyla basıncı artırılır. Basıncı artırılmış hava kütlesi yanma odasında sabit basınç altında yüksek sıcaklık değerlerinde yakılır. Yanma olayı sonucu elde edilmiş gaz kütlesi türbin kademelerinden geçerken genleşir ve türbinin dönme hareketine neden olur. Türbinde meydana Şekil 1: AIRBUS firmasına ait askeri nakliye uçağı A400M’in gelen dönme hareketi aynı zamanda aynı şafta bağlı TP400 Turboprop motor modeli [4] 97 Otomatik Kontrol Ulusal Toplantısı, TOK2013, 26-28 Eylül 2013, Malatya 2. Motor Kontrolü kullanmış olduğu ölçülen değişkenlerden bazıları görülmektedir. Günümüzde değişen ve gelişen teknoloji, hava araçlarını hem yapısal anlamda hem de sistemsel anlamda değişime Tablo 1: Uçak Motorlarında Kullanılan Bazı Kontrol zorlamaktadır. Bu değişimde öncelik ise hız, verim, maliyet Değişkenleri [9] gibi parametrelere dayandırılmaktadır. Uçak motorları da uzun zamandan beri belirtilen bu değişim sürecinde olup, daha önceleri motor kontrol ünitesi dahi olmayan uçaklarda artık gelişmiş bir motor kontrol sistemi kullanılmaktadır. Kontrol Değişken TFE731 F100 F110 GE90 kavramının sisteme dahil olmasıyla birlikte eskiden kullanılan hidromekanik sistemlerin yerini günümüzde elektronik kontrol sistemleri almıştır. FADEC (Full Authority Digital Engine Fan Girişi X X X X Control) adı verilen bu sistemler, uçakta yer alan diğer Toplam sistemlerle koordineli olarak çalışarak durgun durum ve geçiş Sıcaklık hallerinde motor kontrolü sağlamaktadır. Ana fonksiyonu kontrollü ve güvenli uçuş operasyonu için uygun yakıt akışını sağlamaktır [5]. FADEC sistemi [6], hidromekanik sistemlere Fan Girişi X X X göre daha operasyonel esneklik sağlayan ve daha güvenilir bir Basınç sistemdir. sistemlere geçiş yapılmıştır. FADEC sistemi genel olarak iki kısımdan oluşur: Bunlar sırasıyla Elektronik Kontrol Ünitesi (EKÜ) ve Hidromekanik Kompresör X Ünite (HMU)’leridir. HMU kontrol ünitesinden gelen Giriş elektriksel sinyalleri hidromekanik basınca dönüştürerek, Basıncı motor üzerindeki valfleri ve eyleyicileri sürer. Bu iki ünite dışında motor üzerinde pek çok sensör bulunur. Sensörlerin bir kısmı doğrudan kontrol ünitesine bağlı olarak çalışırken, bir Kompresör X diğer kısmı da sistem durumunu görüntülemek için geri Giriş besleme amaçlı olarak kullanılır [7]. Şekil 2’de uçak motoru Sıcaklığı kontrol sistemi basit bir yapıda özetlenmiştir. Kompresör X X X Çıkış Basıncı Kompresör X Bleed Basıncı Egzoz X Çıkış Şekil 2: Uçak Motoru Kontrol Sistemi Sıcaklığı Uçak motoru üreticileri, gaz türbinli motorlardaki kontrol algoritmalarının geliştirilmesi için öncelikli olarak motor Düşük X X X X termodinamik modeline ihtiyaç duymaktadırlar. Bu modeller, Basınç termodinamik yasalar, enerjinin korunumu kanunu ve ampirik Türbini (tecrübeye dayalı) olarak hazırlanan bazı fonksiyonlar Devir kullanılarak oluşturulmakta ve daha sonra hazırlanacak Sayısı kontrol algoritması için temel teşkil etmektedir. Oluşturulan termodinamik modelin verileri kullanılarak sistemde yer alacak eyleyici, dönüştürücü gibi diğer enstrümanların da Yüksek X X X X dahil edilmesiyle hazırlanacak olan tasarım motor Basınç modellemesi ise üç farklı şekilde yapılabilmektedir [8]: Türbini Dinamik Karakteristik, Doğrusal Dinamik Modeller ve Gerçek Devir Zamanlı Modeller’dir. Sayısı Kontrol algoritmaları oluşturulurken kullanılacak kontrol değişkenleri ve ölçülen değişkenler üreticiden üreticiye ve motor tipine göre değişiklik göstermektedir. Kontrol sisteminin komplekslik seviyesi bu parametrelerin ve ölçümü alınan diğer parametrelerin sayısına göre belirlenmektedir [9]. Tablo 1’de dünyanın önde gelen uçak motoru üreticilerinin 98 Otomatik Kontrol Ulusal Toplantısı, TOK2013, 26-28 Eylül 2013, Malatya 3. ANFIS Yapısı x giriş, i düğüm noktası sayısı ve Ai ise bu düğüm noktasının bulanık kümesidir. Her bir düğümde üyelik Doğrusal olmayan sistemlerin modellenmesi, sistemin doğası fonksiyonu olarak en çoğu 1 ve en azı 0 olan çan eğrisi üyelik gereği çok zordur. Bu tip sistemlerin modellenmesinde girişi ve fonksiyonları kullanılır [8]. Fonksiyondaki sabitler ANFIS’de çıkışı haritalayan metotlar kullanmak uygun olmaktadır [10]. öncül (premise) parametreler olarak isimlendirilmektedir [13]. Günümüzde Uyarlanır Sinirsel Bulanık Çıkarım Sistemi 2.Katman: Bu katmanda, düğüm noktasına gelen gelen üyelik (ANFIS) olarak bilinen, bulanık mantık sınıflandırması ile fonksiyonları birbirleriyle çarpılmaktadır. Her düğüm noktası yapay sinir ağları öğrenme mantığının birleştiği yapı giriş- çıktısı bir kural için o kuralın tetikleme veya ateşleme ağırlığı çıkış haritalama işini başarılı bir şekilde yapmaktadır. Bu (firing strength of a rule) olarak tanımlanmaktadır [13]. çalışmada, ANFIS arayüzü yardımıyla pervane devir sayısı, w (x) (y), i1,2. (2) motor devir sayısı ve egzoz çıkış sıcaklığı değerleri baz i Ai Bi alınarak uçak motoru yakıt akış kontrolü modellenmiştir. Daha sonra normalize edilmiş kural tetikleme veya ateşleme Hazırlanan modelde, TEI Tusaş tarafından sağlanan küçük bir ağırlık değerleri ilgili düğüm noktasının tetikleme veya ateşleme turboprop motora ait 3 farklı değişken verisi kullanılmıştır. ağırlık değerinin tüm kuralların tetikleme veya ateşleme ağırlık Kullanılan değişkenler hemen her tip motor için ölçümü alınan değerinin toplamına oranı alınarak denklem (3)’teki gibi kritik değişkenlerdir. Bununla beraber modellemelerde hesaplanır. kullanılacak değişkenler, motor tipine bağlı olarak veya elde w ekstra ölçüm verisi olmasına göre farklılık gösterebilir. wi i , i1,2. (3) w w ANFIS arayüzü temel olarak iki bölümden oluşmaktadır: 1 2 3.Katman: Bu katmanın çıkış fonksiyonu denklem (4)’deki gibi Bunlar sırasıyla FIS (Fuzzy Inherence System) yapısının oluşturulduğu bölüm ve FIS yapısının eğitildiği bölümlerdir. tanımlanmaktadır. {pi, qi, ri} parametreleri ANFIS’de soncul (consequent) parametreler olarak tanımlanmaktadır [13]. FIS yapısının oluşturulduğu bölümde, MATLAB iki çeşit kümeleme metodu sunmaktadır: grid bölümleme (Grid Partitioning) ve eksiltici kümeleme (Subtractive Clustering) Oi3wifi wi(pixqiyri) (4) Ayrıca, daha önceden oluşturulmuş bir FIS yapısı da doğrudan 4.Katman: En son katmanda sistemin çıkışı hesaplanmaktadır. arayüze yüklenerek kullanılabilmektedir. FIS yapısının eğitimi Sistem çıkışı denklem (5)’teki gibi bir önceki düğüm için de yine MATLAB iki yöntem sunmaktadır: Bunlar geri noktasından gelen verilerin tümünün toplanması ile elde yayılım tabanlı öğrenme ve geri yayılım + en küçük kareler edilmektedir. yöntemi kombinasyonundan meydana gelen hibrit öğrenme metodudur [11]. Çalışmamızda kullanmış olduğumuz ANFIS w f i i yapısı aşağıdaki gibidir [10]: Sistem Çıkışı = O4 w f i (5) Modellenecek sistemin x, y şeklinde iki giriş, bir çıkış i i i w i i değeri olduğu varsayılarak ve bu değerler kullanılarak Takagi ve i Sugeno tipi iki kural oluşturulduğunda [12]: Hibrit öğrenme iki prosedürden oluşmaktadır: Bunlar 1. Kural: Eğer x A1 ise ve y B1 ise, çıkış değeri f1=p1x+q1y+r1 sırasıyla İleri Geçiş ve Geri Geçiş Prosedürleri’dir [10]: 2. Kural: Eğer x A2 ise ve y B2 ise, çıkış değeri f2=p2x+q2y+r2 şeklinde olur. Tablo 2: Hibrit Öğrenme Algoritması [10, 11] 1.Katman: Bu katmandaki her düğüm noktasının denklem (1)’deki gibi bir üyelik fonksiyonu mevcuttur. Bu fonksiyon İleri Geçiş Geri Geçiş aynı zamanda düğüm noktası çıkış fonksiyonudur. Öncül Geri yayılım Sabit Oi1A(x) (1) Parametreler Gradyen Azalma i Soncul En Küçük Kareler Sabit Parametreler Tahmini Düğüm Noktası Sinyaller Hata Oranları Çıktıları İleri geçiş prosedüründe, giriş ve ilgili sinyaller 3. katmana kadar bir önceki bölümde bahsedilen işlemlerden geçerek gelir ve en küçük kareler yöntemi aracılığıyla soncul parametreler tahmin edilmeye çalışılır. Geri geçiş prosedüründe ise, gradyen azalma yöntemi kullanılarak çıkıştan girişe tüm düğüm noktalarının hata oranları toplanarak öncül parametrelerin değerleri tekrar güncellenir [10]. 4. Simülasyon Sonuçları Uçak motoru yakıt akışı kontrolü modeli TEI Tusaş Motor Sanayii A. Ş. tarafından tasarlanan küçük bir turboprop motorun motor devir sayısı, pervane devir sayısı, egzoz çıkış sıcaklığı ve yakıt akış verileri kullanılarak Şekil 3: Örnek bir ANFIS Yapısı MATLAB/Simulink üzerinde ANFIS arayüzü aracılığıyla oluşturulmuştur. Elde herhangi bir matematiksel model olmamasına rağmen yalnızca test verileri kullanılarak yakıt 99 Otomatik Kontrol Ulusal Toplantısı, TOK2013, 26-28 Eylül 2013, Malatya akış miktarı hesaplanmaya çalışılmıştır. 981 adet verinin % 80’i eğitim, kalan % 20’lik kısmı ise test verisi olacak şekilde rasgele gruplandırılmıştır. FIS yapısının oluşturulmasında ANFIS arayüzü üzerinde tanımlı eksiltici kümeleme kullanılmıştır. Veri sayısının fazlalığı, artan üyelik fonksiyonu sayısının iterasyon sayısını dolayısıyla işlem süresini uzatması ve artan bu süreye bağlı olarak optimum üyelik fonksiyonu sayısının hesaplanamaması sebebiyle grid bölümleme yöntemi tercih edilmemiştir [14, 15]. Eksiltici kümeleme yöntemi ile en düşük hatayı veren en uygun üyelik fonksiyonu sayısı hesaplanmıştır. Öğrenme metodu olarak da hibrit öğrenme yöntemi kullanılmıştır. Şekil 4’de oluşturulan MATLAB/Simulink modeli, Şekil 5 ve 6’da optimum sayıdaki üyelik fonksiyonuna göre hazırlanmış eğitim verisi ve test verisi orijinal çıktılarının ve ANFIS çıktılarının Şekil 6: Test Verisi ile ANFIS Çıktısının Karşılaştırılması karşılaştırılması, Tablo 3’de ise farklı etki alanı ve üyelik (Etki Alanı: 0.2, Üyelik Fonksiyonu Sayısı:8 için) fonksiyonlarında elde edilen ortalama eğitim ve test veri hataları görülmektedir. Eğitim veri seti için optimum hata her Tablo 3: Eğitim ve Test Verileri Ortalama Hata Tablosu giriş için 10 üyelik fonksiyonu ataması yapıldığında hesaplanmıştır. Test veri seti için optimum hata her giriş için 8 üyelik fonksiyonu ataması yapıldığında hesaplanmıştır. Etki Etki Üyelik Ortalama Ortalama alanı arttıkça üyelik fonksiyonu sayısı azalmaktadır. Bunun Alanı Fonksiyonu Eğitim Verisi Test Verisi nedeni, her verinin çevresinde bulunan veri yoğunluğuna bağlı Sayısı Hatası Hatası olarak potansiyel bir kümelenme merkezi durumunda olması ve etki alanı yarıçapı arttıkça yoğunluğa bağlı olarak kümelenme merkezi seçilen noktanın çevresinde bulunan 0.1 10 4.7722 4.5192 diğer verilerin de bu verinin sınıfına dahil edilecek olması ve böylece toplam üyelik sayısının azalmasıdır [16]. Uygun üyelik fonksiyonu sayısı işleminin bulunmasını ANFIS, etki 0.15 10 4.5011 4.5618 alanı kullanıcı tarafından tanımlandığında otomatik olarak yapmaktadır. 0.2 8 4.7592 4.3686 0.25 7 4.6662 4.394 0.3 7 4.6413 4.3822 0.35 7 4.652 4.3868 0.40 5 5.1591 4.8675 Şekil 4: Simülasyon Modeli 0.45 4 5.8822 5.0731 0.50 3 5.8461 5.0814 0.55 3 5.7024 4.9618 0.60 3 5.5889 4.8737 0.65 3 5.4786 4.7989 0.70 3 5.3952 4.7679 Şekil 5: Eğitim Verisi ile ANFIS Çıktısının Karşılaştırılması (Etki Alanı: 0.15, Üyelik Fonksiyonu Sayısı:10 için) 100 Otomatik Kontrol Ulusal Toplantısı, TOK2013, 26-28 Eylül 2013, Malatya [6] D. Esler, “FADEC's Benefits Today and Tomorrow”, 0.75 3 5.322 4.7404 Business & Commercial Aviation, Cilt: 97, No: 5, 2005. [7] S. Koçak, Gaz Türbinli Motorlarda FADEC Sisteminin Ayrıntılı Olarak İncelenmesi, Lisans Bitirme Tezi, 0.80 3 5.2969 4.7789 Anadolu Üniversitesi, 2012 [8] G.G. Kulikov, H.A. Thompson, Dynamic Modelling of Gas Turbines, Springer, 2005. 0.85 3 5.2744 4.8175 [9] L.C. Jaw, J.D. Mattingly, Aircraft Engine Controls, AIAA Education Series, 2009. [10] J.R. Jang, “ANFIS: Adaptive-Network-Based Fuzzy 0.90 3 5.2549 4.8523 Inference System”, IEEE Transactions on Systems, Man, and Cybernetics, 23, 3, 665-685, 1993. 0.95 3 5.2686 4.9164 [11]http://www.mathworks.com/help/fuzzy/anfis-and-the anfis-editor-gui.html#bq97_i_ [12] T. Takagi, M. Sugeno, “Derivation of fuzzy control rules 1 3 5.2704 4.9632 from human operator’s control actions”, Proceedings of IFAC Symposium on Fuzzy Information, Knowledge Representation and Decision Analysis, 55- 60, 1983. [13] O. Kaynar, M. Zontul, F. Demirkoparan, “Ham petrol Fiyatlarının ANFIS ile Tahmini”, ABMYO Dergisi, 5. Sonuçlar Cilt:17, 3-14, 2010. [14] M.R. Minaz, A. Gün, M. Kurban, N. İmal, “Bilecik İlinin Bu çalışmada, pratikte havacılıkta kullanılan FADEC Farklı Yöntemler Kullanılarak Basınç, Sıcaklık ve Rüzgâr sisteminin sorumlu olduğu yakıt akışı kontrolünün daha Hızı Tahmini”, Gaziosmanpaşa Bilimsel Araştırma sadeleştirilmiş olan yapısı MATLAB/Simulink üzerinde Dergisi, Cilt: 3, 2013. ANFIS arayüzü kullanılarak modellenmiştir. Kontrol parametresi olarak kullanılan giriş verileri egzoz çıkış [15]S. Chiu, “Fuzzy Model Identification Based on Cluster sıcaklığı, pervane devir sayısı ve motor devir sayısı verileri Estimation”, Journal of Intelligent & Fuzzy Systems, 2, 3, kullanılarak uygun yakıt akışı değeri elde edilmeye 1994. çalışılmıştır. Grid bölümleme metodunda üyelik fonksiyonu [16] M. Wei, B. Bai, A. H. Sung, Q. Liu, J. Wang, M. E. sayısı arttıkça iterasyon süresi uzamış ve optimum üyelik Cather, “Predicting injection profiles using ANFIS”, fonksiyonu değeri hesaplanamamıştır. Dolayısıyla diğer Information Sciences, 177, 4445–4461, 2007. yöntem olan eksiltici kümeleme yöntemi farklı etki alanlarına bağlı olarak denenmiş ve farklı sayıda oluşturulan üyelik fonksiyonlarına bağlı olarak optimum üyelik fonksiyonu sayısı hesaplanmıştır. Hesaplamalar sonucunda gerçek motor yakıt akışı değerlerinin oldukça yaklaşık değerlerle ANFIS üzerinde kontrol parametrelerine bağlı olarak elde edilebildiği gözlemlenmiş ve ANFIS üzerinde model için kullanılabilecek en uygun üyelik fonksiyonu değeri tanımlanabilmiştir. Teşekkür Çalışmaya vermiş oldukları teknik destekten dolayı TEI Tusaş Motor Sanayii A. Ş.’ye teşekkür ederiz. Kaynakça [1] H. Karakoç, E.T. Turgut, Gaz Türbinli Motor Sistemleri, T.C. Anadolu Üniversitesi Yayınları, 2008. [2] H. Aydın, Ö. Turan, A. Midilli, T. H. Karakoc, “Transient Performance of an Experimental Turboprop Engine: An Energetic Evaluation”, The Fourth International Conference on Applied Energy (ICAE2012), 2012. [3] Introduction to Turboprop Engine Types http://www.skybrary.aero/bookshelf/books/1626.pdf [4] I. Moir, A. Seabridge, Aircraft Systems, John Wiley & Sons Ltd., 2008. [5] Xicoy Electronica XL, FADEC System-Autostart 06 User Guide. 101 Otomatik Kontrol Ulusal Toplantısı, TOK2013, 26-28 Eylül 2013, Malatya ˘ ˙ ˙ ˙ DOGRUSAL OLMAYAN SISTEMLER IÇIN RUNGE-KUTTA ˙ ˙ MODEL-TABANLI UYARLAMALI PID DENETLEYICI MeriçÇETI˙N1,SelamiBEYHAN2,SerdarI˙PLI˙KÇI˙2 1BilgisayarMühendislig˘iBölümü PamukkaleÜniversitesi,20070,Denizli-Türkiye {mcetin}@pau.edu.tr 2ElektrikElektronikMühendislig˘iBölümü PamukkaleÜniversitesi,20070,Denizli-Türkiye {sbeyhan}@pau.edu.tr,{iplikci}@pau.edu.tr Özetçe tercih edilmektedir. Tasarımı PID denetleyicilerden daha zor Bu çalıs¸mada, Runge-Kutta model-tabanlı model-öngörülü olsadagünümüzdemikrois¸lemcileregömülerekendüstrideçok uyarlamalı Oransal-˙Integral-Türevsel (PID) denetleyici (RK- kullanılmaktadır. Model-öngörülü denetleyicilerin giris¸-çıkıs¸ PID) tanıtılmıs¸ ve dog˘rusal olmayan sistemlerin kontrolünde kestirimufuklarısayesindeileriyeyönelikhataküçüklemesiile kullanılmıs¸tır. Önerilen denetleyici, sürekli-zamanlı bir siste- sistemperformansıoldukçaiyidir.Fakatgradyantemellidenet- min Runge-Kutta ayrıklas¸tırılmıs¸ modeli ile model-öngörülü leyiciolmalarındandolayıgürültüvebozucuetkilerekars¸ıda- performans kriterini azaltacak biçimde PID denetleyici pa- yanıklıdeg˘ildir.MPCyöntemlerininzamandomenindeformule rametrelerini uyarlamaktadır. Sisteme uygulanan giris¸ is¸areti, edilebiliyorolması,ilerleyen-ufuközellig˘ivedurumvekontrol uyarlamalı PID denetleyici ve model-öngörülü denetleyici ile kısıtlarınıkolaycadikkatealabilmegibiüstünözelliklerisayı- eldeedilendüzeltmeterimlerindenolus¸maktadır.ÖnerilenRK- labilir. PID denetleyici, sürekli-zamanlı dog˘rusal-olmayan bioreaktör sistemi üzerinde test edilmis¸tir. Elde edilen denetleme sonuç- 4. mertebeden Runge-Kutta (RK) algoritması, dog˘rulug˘u ları,önerilendenetleyicinindog˘rusal-olmayansistemlerinkont- ve kararlı davranıs¸ları nedeniyle dog˘rusal-olmayan sürekli- rolündeoldukçabas¸arılıoldug˘unugöstermis¸tir. zamanlısistemlerinayrıklas¸tırmasındakullanılır.Runge-Kutta (RK)modeliolarakbilinenbuyapı,uyarlamalıdog˘rusalolma- 1. Giris¸ yan model-öngörülü denetleme, tahminleme, Jacobian hesabı, durum ve parametre kestirimi gibi pek çok uygulamada kul- Kontrolmühendislig˘inde,PIDtüründekidenetleyicilerkolayta- lanılmıs¸tır[9,10,11,12].Buçalıs¸madaisedog˘rusal-olmayan sarımıvegürbüzperformanslarınedeniyledig˘erdenetleyiciler sürekli-zamanlısistemlerindenetlenmesindePIDparametrele- içinde en çok tercih edilen denetleyicidir. Dog˘rusal zamanla rininayarlanmasıiçinsisteminRKmodeliilemodel-öngörülü deg˘is¸meyen sistemlerde PID parametrelerinin tasarımı için li- hataperformansıkullanılarakPIDparametreleriuyarlanmıs¸tır. teratürde pek çok yöntem vardır [1, 2, 3]. Fakat dog˘rusal za- Bununyanındareferansvesistemdinamiklerinindeg˘is¸tig˘inok- manladeg˘is¸ensistemlervereferanssinyalinindeg˘is¸tig˘iuygu- talaraPIDdenetleyicisininyetersizolacag˘ıdüs¸ünülerekmodel- lamalarda,PIDdenetleyiciparametrelerininiyibirreferansta- öngörülüdenetleyiciileüretilmis¸düzeltmeterimitasarlanmıs¸- kibi yapabilmesi için uyarlanması gerekmektedir. Bunun ya- tır. Düzeltme terimi K− adımlık bir tahminleme hatasının en nında dog˘rusal-olmayan sistemler için PID denetleyici tasarı- küçükyapılmasıylaeldeedilir. mındadog˘rullas¸tırmayaparakdengenoktalarındaklasikayar- lama yöntemleri ile parametreler ayarlanmalıdır. Fakat denge Çalıs¸manın devamında: Bölüm 2’de önerilen RK model- noktasısüreklideg˘is¸ebileceg˘indenbuyöntemkullanıs¸lıdeg˘il- tabanlıPIDyapısıdetaylıbirs¸ekildeanlatılacaktır.Tasarlanan dir. dog˘rusal-olmayandenetleyiciyapısınınölçütsistemeuygulan- Dog˘rusallas¸tırmayapmadandog˘rusal-olmayanPIDdenet- masıileeldeedilensonuçlar,Bölüm3’degösterilmektedir.Bö- leme yapılacak ise uyarlamalı PID denetleyicisi kullanılır. Li- lüm4’deiseçalıs¸manınsonuçlarıverilmektedir. teratürde, belirsiz sistemler için kayan-kip uyarlamalı PID [4],matematikselmodelibilinmeyensistemleriçinuyarlamalı yapay-sinirag˘ıileuyarlamalıPID[5]vedestek-vektörmeka- nizmalı uyarlamalı PID [6] denetleyicileri önerilmis¸tir. Bun- 2. ÖnerilenRunge-KuttaModel-Tabanlı ların yanında [7, 8] çalıs¸malarında, PID parametreleri model- UyarlamalıPIDDenetleyici öngörülükontrolmetodukullanılarakamaçfonksiyonununen küçüklenmesiyleayarlanmıs¸tır. Dog˘rusal ve dog˘rusal-olmayan sistemlerin kontrolünde Bubölümde,RKayrıklas¸tırmasıveönerilenRKmodel-tabanlı model-öngörülü denetleyiciler hızlılık ve dog˘ruluk açısından model-öngörülüuyarlanabilirPIDdenetleyicianlatılacaktır. 102 Otomatik Kontrol Ulusal Toplantısı, TOK2013, 26-28 Eylül 2013, Malatya 2.1. Runge-KuttaAyrıklas¸tırması s¸eklindedir.(3)denklemias¸ag˘ıdakigibivektörelformdayazıla- bilir. As¸ag˘ıdakigibikapalıformda ˆx =ˆf(ˆx ,u ), n+1 n n =ˆx +k (5) x˙ =f(x,u), n+1 n ˆy =g(ˆx ,u ). y=g(x,u), (1) n+1 n n burada u∈U,x∈ X,∀t≥0 k +2k +2k +k 11 12 13 14 verilenbirsistemindenetlenebilmesiiçinenuygunkontrolis¸a- 1 k21+2k22+2k23+k24 retininbulunmasıamaçlanmaktadır.Buradakidurumdenklem- kn =6 .. leri . (6) k +2k +2k +k N1 N2 N3 N4 X ={x ∈ℜ|x ≤x ≤x ,i=1,...,N} 1 Ui={u∈i ℜ|u imi≤n u≤iu im}ax (2) = 6(k1+2k2+2k3+k4) min max s¸eklindedir.Buas¸amadansonrat = nT örneklemeanındaki s s¸eklinde kısıtlara bag˘lıdır. Sürekli zamanlı dog˘rusal-olmayan durumdeg˘is¸kenlerix ,··· ,x vegiris¸büyüklüg˘üu veril- 1n Nn n N-boyutlu bir sistemin durum denklemleri vektörel formda dig˘indebirsonrakiörneklemeanındaki(t+T =(n+1)T ) s s (1) s¸eklinde verilmis¸tir. Burada x(t) ∈ ℜN, u(t) ∈ ℜ ve çıkıs¸ı(5)kullanılarakbulunabilir.Artıksürekli-zamanlıbirsis- y(t)∈ ℜolmaküzere,denklemlerdeukontrolis¸aretini,xdu- teminayrık-zamanlıRKyapısındaönerilenuyarlanabilirdenet- rumvektörünü,yçıkıs¸sinyalinigöstermektedir.fivegfonksi- leyicimekanizmasındakullanılabilirhalegelir. yonları;kontrolis¸aretivedurumvektörüdeg˘erlerinegöredina- mikleribilinen,süreklizamanlıvetürevialınabilenfonksiyon- 2.2. Runge-KuttaModel-TabanlıUyarlamalıPID lardır.x ,··· ,x anlıkdurumlarıveu anlıkgiris¸igöster- 1n Nn n meküzeresisteminn.anındakideg˘erleridir.n,t=nT zama- s nındakiörneklemeanınıgösterir.Birsonrakiörneklemeanın- dakisistemeaitx vey s¸eklindeverilendurumveçıkıs¸ in+1 in+1 deg˘erlerias¸ag˘ıdakies¸itliklerdeverilen4.mertebedenRKalgo- ritmasıileöngörülebilir.Bualgoritmadig˘erintegrasyonmetot- larıiçindedahakesinvekararlıoldug˘uiçindog˘rusalolmayan süreklizamanlısisteminayrıklas¸tırılmasıis¸lemindetercihedil- mis¸tir. xˆ =xˆ +k /6+k /3+k /3+k /6, 1n+1 1n 11n 12n 13n 14n . . . xˆ =xˆ +k /6+k /3+k /3+k /6, Nn+1 Nn N1n N2n N3n N4n yˆ =g(xˆ ,...,xˆ ,u ) n+1 1n+1 Nn+1 n (3) burada S¸ekil1:Runge-Kuttamodel-tabanlıPIDdenetleyici k =T f (xˆ ,··· ,xˆ ,u ), 11n s 1 1n Nn n . Önerilen yöntemde, (5) denklemi ile verilen sistemin RK . . modeli, model-tabanlı PID yapısı içinde S¸ekil 1’de görüldüg˘ü kN1n =TsfN(xˆ1n,··· ,xˆNn,un), gibi kullanılmıs¸tır. Burada yˆn modelin çıkıs¸ı, y˜ sistem çıkıs¸ı k =T f (xˆ +0.5k ,··· ,xˆ +0.5k ,u ), y’nin takip etmesi istenen referans is¸areti, en n. örnekleme 12n s 1 1n 11n Nn N1n n anında istenen ve ölçülen çıkıs¸ arasındaki hata, u kontrol n+1 . .. is¸aretiveδun+1 düzeltmeterimidir.UyarlamalıRK-PIDdön- güsü, RK modeli, uyarlamalı PID blog˘u ve iki adet Jacobian k =T f (xˆ +0.5k ,··· ,xˆ +0.5k ,u ), N2n s N 1n 11n Nn N1n n blog˘undan olus¸ur. Yapı içindeki Jacobian bloklarından birisi k13n =Tsf1(xˆ1n +0.5k12n,··· ,xˆNn +0.5kN2n,un), kontrol sinyalini düzeltir, dig˘eri ise parametre ayarlamasında .. kullanılır.PIDdenetleyici,en hatasınıbarındıranveas¸ag˘ıdaki . formdaverilenu s¸eklindekikontrolis¸aretiniüretir. n+1 k =T f (xˆ +0.5k ,··· ,xˆ +0.5k ,u ), N3n s N 1n 12n Nn N2n n un+1 =un+KP(en−en(cid:0)1)+KIen k =T f (xˆ +k ,··· ,xˆ +k ,u ), (7) 14n s 1 1n 13n Nn N3n n +KD(en−2en(cid:0)1+en(cid:0)2) . . . Burada K , K ve K en uygun biçimde ayarlanması ge- P I D k =T f (xˆ +k ,··· ,xˆ +k ,u ). rekenPIDparametreleridir.Buparametrelerindenetlemeis¸le- N4n s N 1n 13n Nn N3n n (4) minebas¸lamadanöncekibas¸langıçkos¸ullarısıfırdır.PIDpara- 103 Otomatik Kontrol Ulusal Toplantısı, TOK2013, 26-28 Eylül 2013, Malatya metreleribas¸langıçdeg˘erleriiledenetlemeis¸lemigerçekles¸tiri- yetersizkalabilir.Bazımodellemehatalarıveharicibozucula- lemedig˘iiçinenuygundeg˘erlerineayarlanmalıdır.Buamaçla, rın kontrol hareketinin üzerindeki olumsuz etkilerden kurtul- sistemin RK modeli kullanılır. Bu model, sistemin K− adım makiçinδu gibibirdüzeltmeterimineihtiyaçduyulur.RK n+1 sonraki gelecek davranıs¸larını tahmin ettig˘i biles¸enlerden olu- model-tabanlıPIDyapısıiçindekidüzeltmeblog˘ununamacı(8) s¸an[yˆ ,yˆ ,...,yˆ ]s¸eklindebiryörüngevektörüüre- denklemindeki amaç fonksiyonunu en küçük yapacak uygun n+1 n+2 n+K tir.Buvektörolus¸turulurkenPIDdenetleyicisiileüretilenu δu terimini bulmaktır. Farklı bir deyis¸le, düzeltme blog˘u n+1 n+1 kontrol is¸areti RK modeline K kere uygulanır. Bu tahminlere amaçfonksiyonuF’yi2.derecedenTayloryaklas¸ımınadaya- dayanarakPIDdenetleyisininparametreleri,K−adımsonraki narakδu teriminegöreenküçükyapmayaçalıs¸ır: n+1 tahminleme hatalarının kareleri toplamını en küçük yapmayı amaçlayan bir maliyet fonksiyonu ile belirlenir. Dig˘er bir de- F(un+1+δun+1)∼=F(un+1)+ ∂∂F(u(un+1))δun+1 yis¸le as¸ag˘ıdaki formda verilen amaç fonksiyonu F en küçük n+1 (14) yapılmayaçalıs¸ılır. + 1∂2F(un+1)(δu )2 2 ∂u2 n+1 ∑K n+1 1 1 F(un+1)= 2 (y˜n+k−yˆn+k)2+ 2λ(un+1−un)2 Amaçfonksiyonunuenküçükyapanδun+1terimiarandıg˘ıiçin k=1 F’ninδun+1’egöretürevialınıpsıfıraes¸itlenirse; (8) burada K kestirim ufku ve λ cezalandırma terimidir. PID pa- @F rametreleriningüncellenmesiveamaçfonksiyonununenküçük δun+1 =−@@u2nF+1 (15) yapılmasıiçinLevenberg-Marquardt(LM)kuralıkullanılmıs¸tır. @u2 n+1 ∆KPID =−(JTJ+µI)(cid:0)1JTeˆ (9) (15)denkleminegöreamaçfonksiyonununδun+1’egöre1.ve 2. türevlerinin hesaplanması gerekmektedir. Ancak 2. derece- K K dentürevlerinhesaplamakarmas¸ıklıg˘ınedeniylebutürevlerye- P P K ⇐= K +∆K (10) rineJacobianyaklas¸ıklıklarıkullanılabilir.(K+1)×1boyutlu I I PID KD KD JmJacobianmatrisi, briutmraadlaarµı apraarsaımndeatreuszilaSs¸tmeeapessatg˘-DlaeysacnenbtirveteGrimaudsisr-.NIe,w3to×n3algboo-- Jm =−[ @@uyˆnn++11 @@uyˆnn++21 ··· @@yˆunn++K1 √λ ]T yutundaki birim matris olup J ise (K +1)×3 boyutundaki (16) Jacobianmatrisidir. s¸eklindeise1.ve2.derecedentürevlerJmvektörükullanılarak as¸ag˘ıdakis¸ekildetanımlanır. @eˆn+1 @eˆn+1 @eˆn+1 @@eˆKn+P2 @@eˆKn+I2 @@eˆKn+D2 ∂F∂u(unn++11) =2JTmeˆ ve ∂2∂Fu(2nu+n1+1) ∼=2JTmJm (17) J= @K...P @K...I @K...D Buas¸amadansonradüzeltmeterimias¸ag˘ıdakigibiyazılabilir. @eˆn+K @eˆn+K @eˆn+K JTeˆ @KP @KI @KD δu =− m (18) n+1 JTJ p p p m m @ (cid:21)(un+1(cid:0)un) @ (cid:21)(un+1(cid:0)un) @ (cid:21)(un+1(cid:0)un) @KP @KI @KD (18)ifadesindendegörüldüg˘ügibiartık1.derecedentürevlere (11) ihtiyaçduyulur.Budurumda(12)denklemindekiJacobianmat- @yˆn+1 @yˆn+1 @yˆn+1 @KP @KI @KD risizincirkuralındanfaydalanılarakas¸ag˘ıdakigibi2farklımat- riss¸eklindeyazılabilir. =− @@@@yˆyˆnKKn...++PPK2 @@y@@ˆyˆnKKn...++IIK2 @@@@yˆyˆKKnn...++DDK2 @@@@uuyyˆˆnnnn++++2111 [ ] veeˆtahmin√haλta@l∆a@rKuınPv+e1ktör√üoλlm@∆a@ukKnüI+z1ere:√λ@∆@KunD+1 (12) J=− @@yˆunn...++K1 @@uKn+P1 @@uKn+I1 @@uKn+D1 (19) √ eˆ y˜ −yˆ λ n+1 n+1 n+1 .. .. =JmJc eˆ= . = . (13) √eˆn+K √y˜n+K−yˆn+K burada Jc vektörü, un+1’in PID parametrelerine göre parçalı λ∆u λ(u −u ) türevleridir.Butürevlersadeceizleyicihatasıkullanılarakas¸a- n+1 n+1 n g˘ıdakiformdayazılabilir: Bu as¸amada artık PID parametreleri her adımda güncellenir. Güncellenenparametrelerherzamankabuledilebilirbirkont- en−en(cid:0)1 T rol üretebilecek kadar uygun bir deg˘erde olmayabilir. Bu du- J = e (20) c n rumdabukontrolhareketisistemiistenenyörüngeyegötürmede en−2en(cid:0)1+en(cid:0)2 104 Otomatik Kontrol Ulusal Toplantısı, TOK2013, 26-28 Eylül 2013, Malatya Sonuçta, J vektörü, PID parametrelerini güncellemek için s¸eklindedir.Ayrıca(24)denklemiiçingerekenkontrolis¸aretine m (12)denklemindekiJacobianmatrisininhesabındaveayrıcadü- göreRKmodel-tabanlıparçalıtürevleriseas¸ag˘ıdakigibidir: zeltme terimini bulmada (15) kullanılır. Uygun bir δu he- n+1 (cid:12) dsaepnladnodlıag˘yıındJa amveakçtöforünnküsniyoynaupnıuiçeinndkeükçiüektkyisaipmaçaıköçzaelglöig˘rüinl-- @u@nk+11 =Ts@@fu(xn;+u1) =Ts@u@nf+1(cid:12)(cid:12)x=x (29) m mektedir.Buçalıs¸mada,J vektörüçalıs¸masındaönerilenRK modeliilehesaplanmaktadmır. @u@nk+21 =Ts@f(x@+u0n:+51k1;u) RKmSiostdeemliisnaKye−sinaddeımdeinlekrlieymed(5o)g˘irlueoitlearnaçtiıfkbıs¸ilras¸reınkıinldkeegsetirrçimeki-, =Ts(@f@((xx++00:5:5kk11;)u)@(x@+un0+:51k1) + @u@nf+1(cid:12)(cid:12)(cid:12)x=x+0:5k1) ( (cid:12) (cid:12) ) lvees¸ktitröirlüirn.üBnu(raudna)(d[etg˘+is¸mTsetdi+g˘iKvaTrssa]y)ıplmeraikytoadduırn.dBaöaydleacyekiolenrtiryoel =Ts 0:5@@xf(cid:12)(cid:12)x=x+0:5k1@u@nk+11 + @u@nf+1(cid:12)(cid:12)x=x+0:5k1 (30) yönelikkestirimler[yˆ ,...,yˆ ]s¸eklindeolur.J Jaco- bkiuarnalıvyekartödrımün∂ıüyyˆlaolgues¸ntuenrl+aon∂1laTr@@gauyˆk(cid:12)(cid:12)nn++11n,+.K..,@∂@yˆuxˆnn++K1 terimlermi zincir @u@nk+31 ==TTss((@f@(f@x((x@+xu++(cid:12)0n:+050:1k5:52kk;22u;))u)@(x@+un0+:51k2) + @(cid:12)u@nf+1(cid:12)(cid:12)(cid:12)x=x+)0:5k2) ∂unn++K1 = ∂x (cid:12)x=xˆ[n+k]∂unn++k1, (21) =Ts 0:5@@xf(cid:12)(cid:12)x=x+0:5k2@u@nk+21 + @u@nf+1(cid:12)(cid:12)x=x+0:5k2 (31) es¸esakslinhdeseahpelasnapmlaansıırg.e(r2e1k)eens¸i@tlxˆign˘+inkdeifa@@dTxegletreiraims¸aig˘sıdabakititigr.ibBiuardaıdma @u@nk+41 =Ts(@f(x@u+n+k13;u) (cid:12) ) adımbulunur. @un+1 =Ts @f@((xx++kk33;)u)@(@xu+n+k13) + @u@nf+1(cid:12)(cid:12)x=x+k3 (32) ( (cid:12) (cid:12) ) ∂xˆn+k = ∂^f(cid:12)(cid:12)(cid:12) ∂xˆn+k(cid:0)1 + ∂^f (cid:12)(cid:12)(cid:12) =Ts @@xf(cid:12)(cid:12)x=x+k3@u@nk+31 + @u@nf+1(cid:12)(cid:12)x=x+k3 ∂u ∂x x=xˆ[n+k(cid:0)1] ∂u ∂u x=xˆ[n+k(cid:0)1] n+1 un+1=un+1 n+1 n+1 un+1=un+1 s¸eklindedir. (22) burada 3. ÖlçütSistemveBenzetimSonuçları @^f @(x+k) = RK-PID yapısı, MATLAB benzetim ortamında dog˘rusal- @xn @xn (23) =I+ 1@k1 + 1@k2 + 1@k3 + 1@k4: olmayan, sürekli-zamanlı ve matematiksel dinamikleri bilinen 6@xn 3@xn 3@xn 6@xn bioreaktörsistemineuygulanmıs¸tır.Dog˘rusal-olmayanbioreak- törsistemdinamikleri, ve @^f @(x+k) @k c˙ (t)=−c (t)u(t)+c (t)(1−c (t))ec2(t)=(cid:13), = = 1 1 1 2 @un+1 = 16@@uu@nnk++111 + 13@@uu@n+nk+211 + 13@u@nk+31 + 16@u@nk+41 (24) c˙2(t)=−c2(t)u(t)+c1(t)(1−c2(t))ec2(t)=(cid:13)1+1β+−βc(23(3t)) s¸eklindedir.(23)denklemindekidurumlaragörehesaplananRK s¸eklindedir.Buradau(t)akıs¸debisinigöstermeküzerekontrol model-tabanlıparçalıtürevler: is¸aretidir. c1(t) hücre yog˘unlug˘u ve c2(t) birim hacim bas¸ına düs¸enbesinmiktarıolupy(t) = c (t)s¸eklindedir.Bioreaktör (cid:12) 1 @@xkn1 =Ts@f@(xxn;u) =Ts@@xfn(cid:12)(cid:12)uxnn+=1x=nun+1 (25) sdiisrt.eKmoinntirnolnois¸marientainlipnatreapmeedtreeg˘leerril,erγiu=mi0n.4=8,0βve=u0m.a0x2=gib2i-, örneklemesüresiiseτ = τ = 0.01saralıg˘ındatutul- @@xkn2 =Ts@f(x+@x0:n5k1;u) mus¸tur. RK-PID denetmleiynicisinimnaaxmacı, akıs¸ debisini kontrol =Ts@(f(@xx++00::55kk11;u)@(x+@x0n:5k1) ) (26) ederÇekalrıes¸fmeraadnasbies¸anrzeettiinmitsaoknipuçeltamrıebkitrikr.açkestirimufkuiçinde- =Ts @@xfn(cid:12)(cid:12)(cid:12)uxnn+=1x=nu+n0+:51k1(I+0:5@@xkn1) npeanrammies¸ttirrelaenrciaikçienneilydiesoenduilçmKis¸ti=r. K10v<ec1e0zadteeg˘reimrleirλi i=çin0.i0le0r1i yönlü düzeltmenin kısa olmasından dolayı PID parametreleri @@xkn3 =Ts@f(x+@x0:n5k2;u) dog˘ru bir s¸ekilde ayarlanamamaktadır. Bu durum çıkıs¸ is¸are- =Ts@(f(@xx++00::55kk22;u)@(x+@x0n:5k2) ) (27) tindeKsa>lını1m0ladreag˘enreldeerinioçilnmiaykitlaeds¸tıirr.meolmaklaberaberönemli =Ts @@xfn(cid:12)(cid:12)(cid:12)uxnn+=1x=nu+n0+:51k2(I+0:5@@xkn2) olarnaancdaakbtiürrdeevg˘iiss¸¸ilmemyloerkituarr.ttBıg˘uıniuçninyhanesınadpalaKma−saüdreımsiiuçzinamheaskatpa-- dır.BuyüzdenenuygunKdeg˘eri10olarakseçilmis¸tir.λpara- @@xkn4 =Ts@f(x@+xnk3;u) metresinin0.001’denküçükdeg˘erleriiçinkontrolis¸aretindeki =Ts@(f(@xx++kk33;u)@(x@+xnk3) ) (28) cteedziarl.aλndı>rma0.k0ü0ç1ükdeog˘lderulg˘euriiiçlienddaahhaagheızçlıizilzelmemeeyaypaıplaılbmilamseınka- (cid:12) rag˘men kontrol is¸aretinde osilasyona neden oldug˘undan siste- =Ts @@xfn(cid:12)(cid:12)uxnn+=1x=nu+nk+31(I+ @@xkn3) mindahageçoturmasınanedenolmaktadır. 105
Description: