ΠΤΥΧΙΑΚΗ ΕΡΓΑΣΙΑ ΘΕΜΑ:

1. ΤΕΧΝΟΛΟΓΙΚΟ ΕΚΠΑΙΔΕΥΤΙΚΟ ΙΔΡΥΜΑ ΣΕΡΡΩΝΣΧΟΛΗ ΤΕΧΝΟΛΟΓΙΚΩΝ ΕΦΑΡΜΟΓΩΝΤΜΗΜΑ ΠΛΗΡΟΦΟΡΙΚΗΣ & ΕΠΙΚΟΙΝΩΝΙΩΝ ΠΤΥΧΙΑΚΗ ΕΡΓΑΣΙΑ ΘΕΜΑ: Ανάλυση και σύνθεση ενός ανάστροφου εκκρεμούς χρησιμοποιώντας τεχνικές ευφυή αυτομάτου ελέγχου με χρήση του MATLAB Σπουδάστρια: Χρυσούλα Γκέγκα

2. ΠΕΡΙΓΡΑΦΗ ΤΟΥ ΠΡΑΓΜΑΤΙΚΟΥ ΣΥΣΤΗΜΑΤΟΣ ΑΝΑΣΤΡΟΦΟΥ ΕΚΚΡΕΜΟΥΣ • Τα μέρη από τα οποία αποτελείται το πραγματικό μοντέλο του ανάστροφου εκκρεμούς είναι: α) ένα βαγόνι τοποθετημένο πάνω σε δύο ράγες και β) μία ράβδο σε όρθια θέση η οποία είναι βιδωμένη πάνω στο βαγόνι και μπορεί να κάνει τις ίδιες κινήσεις με αυτές του βαγονιού

3. ΓΙΑΤΙ ΜΑΣ ΕΝΔΙΑΦΕΡΕΙ Η ΜΕΛΕΤΗ ΤΟΥ • Μας ενδιαφέρει γιατί μπορεί να βρει εφαρμογή σε αρκετά συστήματα που έχουν παρόμοια λειτουργία με το εκκρεμές, μερικά από αυτά είναι: • Το αρχικό στάδιο πτήσης ενός πύραυλου. • Ο ρομποτικός βραχίονας των robot’s • Η ισορροπία ακόμη και του ίδιου του ανθρώπου τόσο όταν είναι σε κίνηση όσο και όταν παραμένει στάσιμος σε ένα σημείο.

4. ΣΤΟΧΟΙ ΤΗΣ ΕΡΓΑΣΙΑΣ • Το όλο σύστημα να έρχεται στην επιθυμητή θέση ανάλογα με την είσοδο που του δίνεται κάθε φορά. • Η ράβδος θα πρέπει να παραμένει όρθια σε όποια ενέργεια και αν εκτελείται από το σύστημα κάθε φορά. • Το σύστημα θα πρέπει να αντεπεξέρχεται κάτω από διάφορους άλλους «ανεπιθύμητους» παράγοντες που δεν ήταν αναμενόμενοι όπως ο θόρυβος και οι διάφορες αλλαγές των παραμέτρων.

5. TO ΜΑΘΗΜΑΤΙΚΟ ΜΟΝΤΕΛΟ ΤΟΥ ΑΝΑΣΤΡΟΦΟΥ ΕΚΚΡΕΜΟΥΣ • Το ανάστροφο εκκρεμές αποτελείται από ένα βαγόνι που κινείται πάνω σε μία ράγα, ενωμένο με μία κινητή ράβδο πουπεριστρέφεται γύρω από έναν άξονα • Το διάνυσμα x= col(x1, x2,x3, x4) περιγράφει την κατάσταση του συστήματος, όπου x1: η θέση του βαγονιού, x2: η γωνία του εκκρεμούς, x3: η ταχύτητα του βαγονιού και x4: η γωνιακή ταχύτητα του εκκρεμούς. • Στο βαγόνι εφαρμόζεται μία ελεγχόμενη δύναμη F -παράλληλη στη ράγα-, η οποία ελέγχεται από ένα μοτέρ ηλεκτρικής τάσης u, όπου το u κυμαίνεται στο [-0.5, 0.5]

6. ΟΙ ΕΞΙΣΩΣΕΙΣ ΤΟΥ ΣΥΣΤΗΜΑΤΟΣ

7. mpw : μάζα του φορτίου [kg] mps : μάζα του πόλου [kg] lp: μήκος του πόλου [m] lpo: απόσταση μεταξύ του κέντρου της μάζας του πόλου και του άξονα περιστροφής [m] lpwo: απόσταση μεταξύ του κέντρου της μάζας του φορτίου και της περιστροφής του άξονα του εκκρεμούς [m]. lc: μήκος του φορτίου [m]. lco: απόσταση μεταξύ του κέντρου της μάζας φορτίων και της περιστροφής του άξονα του εκκρεμούς [m]. rp: ακτίνα του πόλου [m]. rc: ακτίνα του φορτίου [m] mc: ισοδύναμη μάζα του βαγονιού, των τροχαλιών και του συνεχούς στροφέα. m: ισοδύναμη μάζα εκκρεμούς και κάρου l: απόσταση από τον άξονα περιστροφής και το κέντρο της μάζας του συστήματος fc: δυναμικός συντελεστής τριβής του κάρου fs: στατικός συντελεστής τριβής του κάρου fp:περιστροφικός συντελεστής τριβής Jp: στιγμή αδράνειας του εκκρεμούς όσον αφορά τον άξονα περιστροφής g: βαρύτητα p1: δύναμη ελέγχου στην αναλογία σημάτων PWM p2: δύναμη ελέγχου στην αναλογία ταχύτητας του κάρου umax: μέγιστη παραγόμενη τιμή σήματος από το PWM R1: μήκος της ράγας T: περίοδος του εκκρεμούς J: στιγμή της αδράνειας OI ΠΑΡΑΜΕΤΡΟΙ ΤΟΥ ΣΥΣΤΗΜΑΤΟΣ

8. TOSIMULINK ΜΟΝΤΕΛΟ

9. ΕΥΦΥΗΣ ΕΛΕΓΧΟΣ • Στηρίζεται στην ανθρώπινη γνώση και εμπειρία • Η γνώση αναπαρίσταται από λεκτικούς κανόνες (Linguist Rules) που είναι της μορφής “εάν (αίτια) τότε (συμπέρασμα)” “If x then y” • Μέρη του αποτελούν η Ασαφής Λογική, τα Τεχνητά Νευρωνικά Δίκτυα, τα Έμπειρα Συστήματα καθώς και οι συνδυασμοί τους (Υβριδικά Συστήματα) • Χρησιμοποιείται σε πολλές βιομηχανικές εφαρμογές όπως σε βιομηχανίες τσιμέντου, πετρελαιοειδών, λιπασμάτων

10. ΠΛΕΟΝΕΚΤΗΜΑΤΑ-ΜΕΙΟΝΕΚΤΗΜΑΤΑ ΕΥΦΥΗ ΕΛΕΓΧΟΥ ΠΛΕΟΝΕΚΤΗΜΑΤΑ • Στηρίζεται στη μίμηση του ανθρώπου και την αναπαραγωγή της ανθρώπινης γνώσης • Χρησιμοποιείται σε πολύπλοκα συστήματα στα οποία είναι δύσκολο ή και αδύνατο να βρεθεί το μαθηματικό μοντέλο • Οι ευφυείς ελεγκτές έχουν την ικανότητα να λειτουργούν κάτω από ένα περιβάλλον ασάφειας και αβεβαιότητας ΜΕΙΟΝΕΚΤΗΜΑΤΑ • Οι ευφυείς ελεγκτές δεν κατέχουν την ικανότητα προσαρμογής και μάθησης νέων κανόνων

11. ΑΣΑΦΗΣ ΛΟΓΙΚΗ (FUZZY LOGIC) • Στηρίζεται στη θεωρία των ασαφών συνόλων(θεμελιωτής L. A. Zadeh, 1965) • Χρησιμοποιεί την έννοια του βαθμού συµµετοχής (degree of membership) • Αναπαριστά ανακριβείς ή αβέβαιες γνώσεις (π.χ. ο Γιώργος είναι ψηλός) • Εφαρμόζεται στην παραγωγή φωτογραφικών μηχανών, βιντεοκαμερών, κλιματιστικών, πλυντηρίων κτλ.

12. ΑΣΑΦΗ ΣΥΝΟΛΑ • Ένα ασαφές σύνολο (fuzzy set) , έστω Α, ορίζεται µέσω της συνάρτησης συμμετοχής (membership function): Το Χ είναι το υπερσύνολο αναφοράς (universe of discourse) και x ένα στοιχείο του Χ. • Το µΑ(x) είναι ένας πραγματικός αριθμός (0 ≤µΑ(x) ≤1) που δείχνει τον βαθμό κατά τον οποίο το x ανήκει στο Α και ονομάζεται βαθμός συµµετοχής(degree of membership) • Σε περίπτωση συνεχών τιμών, το Α εκφράζεται σαν μια συνάρτηση όπως triangular, trapezoidal, gaussianκτλ.

13. ΠΛΕΟΝΕΚΤΗΜΑΤΑ-ΜΕΙΟΝΕΚΤΗΜΑΤΑ ΤΗΣ ΑΣΑΦΗΣ ΛΟΓΙΚΗΣ ΠΛΕΟΝΕΚΤΗΜΑΤΑ • Ανεκτικότητα σε ανακριβή γνώση • Ευκολία στην αναπαράσταση της γνώσης με τη χρήση λεκτικών όρων (π.χ. κοντός, μέτριος, ψηλός) • Ευκολία στην κατανόηση και στη χρήση καθώς βασίζεται στη φυσική γλώσσα ΜΕΙΟΝΕΚΤΗΜΑΤΑ • Δυσκολία στο σχεδιασμό ενός ασαφή ελεγκτή (fuzzy controller)εφόσον δεν υπάρχει κάποια συγκεκριμένη θεωρία • Απαραίτητη η γνώση ενός εμπειρογνώμονα

14. ΑΣΑΦΕΣ ΜΟΝΤΕΛΟ ΤΥΠΟΥ MAMDANI 1.Ασαφοποίηση(Fuzzification) Μετατροπή των πραγματικών εισόδων σε ασαφή σύνολα 2. Εκτίμηση κανόνων (Rule evaluation) 2.1 Εκτίμηση συνθηκών Υπολογισμός του συνδυασμένου βαθμού συµµετοχής 2.2 Εκτίμηση συµπεράσµατος Εφαρμογή του συνδυασμένου βαθμού στη συνάρτηση συµµετοχής του συµπεράσµατος 3. Συνάθροιση εξόδων (Aggregation) Συνδυασμός των παραχθέντων συναρτήσεων συµµετοχής των εµπλεκόµενωνκανόνων σ’ ένα ασαφές σύνολο 4. Αποασαφοποίηση (Defuzzification) Μετατροπή του αποτελέσματος της συνάθροισης σε πραγματική τιμή

15. ΤO SIMULINK ΜΟΝΤΕΛΟ ΤΟΥ ΑΝΑΣΤΡΟΦΟΥ ΕΚΚΡΕΜΟΥΣ ΜΕ ΤΟΝ ΑΣΑΦΗ ΕΛΕΓΚΤΗ

16. ΕΦΑΡΜΟΓΗ ΤΗΣ ΜΕΘΟΔΟΥ MAMDANI ΣΤΟΝ ΑΣΑΦΗ ΕΛΕΓΚΤΗ Είσοδοι του ελεγκτή: α) Cart Position [-2 2] β) Angle [-0.3 0.3] γ) Cart Velocity[-2 2] δ) Angular Velocity[-2 2] Έξοδος του ελεγκτή: Force [- 3 3]

17. ΕΦΑΡΜΟΓΗ ΤΗΣ «ΓΩΝΙΑΣ» ΩΣ ΕΙΣΟΔΟ ΣΤΟΝ ΕΛΕΓΚΤΗ Πραγματοποιήθηκαν πολλές δοκιμές για να βρεθούν οι συνθήκες εκείνες που θα ισορροπούσαν καλύτερα το σύστημα. Οι συνθήκες αυτές που βρέθηκαν είναι: AND (prod) OR (probor) Συνεπαγωγή prod Συνάθροιση probor Αποσαφοποίηση centroid Ζυγός αριθμός ασαφών συνόλων Η μορφή ασαφών συνόλων gauss2mf Επόμενο βήμα ήταν η εύρεση του αριθμού για την κάθε μία από τις εισόδους και για την έξοδο όπωςκαι περισσότεροι κανόνες που θα ισορροπούσαν το σύστημα Κανόνες: If (angle is negative) then (force is very negative) If (angle is positive) then (force is very positive)

18. ΕΦΑΡΜΟΓΗ ΤΗΣ «ΓΩΝΙΑΣ» ΚΑΙ ΤΗΣ «ΓΩΝΙΑΚΗΣ ΤΑΧΥΤΗΤΑΣ» Εφαρμόζοντας τις προηγούμενες συνθήκες και πραγματοποιώντας πολλές αλλαγές στο σύνολο των ασαφών συνόλων των εισόδων, το σύστημα παρατηρήθηκε να είναι πιο ευσταθές ενώ δόθηκαν και περισσότεροι κανόνες Μονός αριθμός ασαφών συνόλων μορφής gauss2mf Ζυγός αριθμός ασαφών συνόλων μορφής gauss2mf • Κανόνες: • If (angle is negative) and (angular velocity is negative) then (force is very negative) • If (angle is negative) and (angular velocity is zero) then (force is negative) • If (angle is negative) and (angular velocity is positive) then (force is small negative) • If (angle is zero) and (angular velocity is negative) then (force is negative) • If (angle is zero) and (angular velocity is positive) then (force is positive) • If (angle is positive) and (angular velocity is negative) then (force is small positive) • If (angle is positive) and(angular velocity is zero) then (force is positive) • If (angle is positive) and (angular velocity is positive) then (force is very positive)

19. ΕΦΑΡΜΟΓΗ ΟΛΩΝ ΤΩΝ ΕΙΣΟΔΩΝ Δοκιμάζοντας εκ νέου αλλαγές το σύνολο των ασαφών συνόλων των εισόδων καθώς και της εξόδου και με την εφαρμογή νέων κανόνων, το σύστημα ισορρόπησε ικανοποιητικά. Ζυγός αριθμός ασαφών συνόλων για κάθε είσοδο μορφής gauss2mf Μονός αριθμός ασαφών συνόλων μορφής gauss2mf

20. ΚΑΝΟΝΕΣ • Οι κανόνες που συντάχθηκαν ήταν της μορφής “If x is A then y is B” • Ο τελεστής OR χρησιμοποιείται για να δηλώσει είτε το μέγιστο (max)είτε το αλγεβρικό άθροισμα (probor)των ασαφών συνόλων Α και Β • Ο τελεστής AND χρησιμοποιείται για να δηλώσει είτε το ελάχιστο (min) είτε το γινόμενο (product) των ασαφών συνόλων Α και Β • Οι κανόνες δημιουργήθηκαν από την γνώση και την προσωπική εμπειρία του σχεδιαστή του συστήματος • Το σύνολο των κανόνων ήταν αποτέλεσμα του συνόλου των εισόδων υψωμένο στον αριθμό των ασαφών συνόλων των εισόδων του ασαφή ελεγκτή. Συγκεκριμένα ο αριθμός των ασαφών συνόλων για την κάθε μία από τις 4 εισόδους του ασαφή ελεγκτή ήταν 2, οπότε 24=16 κανόνες

21. ΚΑΝΟΝΕΣ

22. ΚΑΝΟΝΕΣ 1ος κανόνας If (cart position is negative) AND (angle is negative) AND (cart velocity is negative) AND (angular velocity is negative) THEN (force is very negative) 5ος κανόνας If (cart position is negative) AND (angle is positive) AND (cart velocity is negative) AND (angular velocity is negative) THEN (force is small negative)

23. ΕΦΑΡΜΟΓΗ ΚΑΛΥΤΕΡΟΥ ΕΛΕΓΚΤΗ(1) Είσοδος «Cart Position» Είσοδος «Angle» Έξοδος «Force» 2 ασαφή σύνολα (negative, positive) Είσοδος«Cart Velocity» Είσοδος«Angular Velocity» 7 ασαφή σύνολα (very negative, negative, small negative, zero, small positive, positive, very positive)

24. ΕΦΑΡΜΟΓΗ ΚΑΛΥΤΕΡΟΥ ΕΛΕΓΚΤΗ(2) Είσοδοι Ελεγκτή Έξοδος Ελεγκτή Μέθοδος Συνεπαγωγής Γινομένου Ο τελεστής ANDως γινόμενο Μέθοδος Συνάθροισης Αλγεβρικού Αθροίσματος Ο τελεστής OR ωςαλγεβρικό άθροισμα Μέθοδος Αποσαφοποίησης Κεντρώου

25. ΕΦΑΡΜΟΓΗ ΚΑΛΥΤΕΡΟΥ ΕΛΕΓΚΤΗ(3) Κανόνες Cart Position – Cart Velocity N: negative P: positive Z: zero SN: small negative SP: small positive VN: very negative VP: very positive Angle – Angular Velocity If (cart position is negative) AND (angle is negative) AND (cart velocity is positive) AND (angular velocity is positive) THEN (force is small negative)

26. ΕΦΑΡΜΟΓΗ ΚΑΛΥΤΕΡΟΥ ΕΛΕΓΚΤΗ(4)

27. ΕΦΑΡΜΟΓΗ ΚΑΛΥΤΕΡΟΥ ΕΛΕΓΚΤΗ(5) Πραγματικήθέση του βαγονιού Επιθυμητή θέση του βαγονιού

28. ΕΙΣΑΓΩΓΗ ΘΟΡΥΒΟΥ ΣΤΟ ΣΥΣΤΗΜΑ • Εισήχθη θόρυβος στο σύστημα για να παρουσιαστεί η ευρωστία του ασαφή ελεγκτή.

29. ΕΙΣΑΓΩΓΗ ΘΟΡΥΒΟΥ ΣΤΟ ΣΥΣΤΗΜΑ(2) Ο θόρυβος που εισήχθη σε αυτή τη δοκιμή ήταν πολύ μεγάλος, με αποτέλεσμα ο ελεγκτής να μην μπορέσει να ανταπεξέρθει • Οι τιμές θορύβου στις εισόδους ήταν: α) «Cart Position» ίση με 35·10-3 β «Angle»ίση με 6 ·10-4 γ) «Cart Velocity»ίση με 26 ·10-3 δ) «Angular Velocity»ίση με 8 ·10-2

30. ΕΙΣΑΓΩΓΗ ΘΟΡΥΒΟΥ ΣΤΟ ΣΥΣΤΗΜΑ(3) • Παρουσία θορύβου μειωμένη σε σχέση με την προηγούμενη δοκιμή Οι τιμές θορύβου για τις εισόδους ήταν: α) «Cart Position»ίση με 10-4 β) «Angle» ίση με 10-5 γ) «Cart Velocity» ίση με 10-4 δ) «Angular Velocity» ίση με 10-3

31. ΣΥΜΠΕΡΑΣΜΑΤΑ ΤΗΣ ΕΠΙΡΡΟΗΣ ΤΟΥ ΘΟΡΥΒΟΥ ΣΤΗΝ ΙΣΟΡΡΟΠΙΑ ΤΟΥ ΣΥΣΤΗΜΑΤΟΣ • Το σύστημα του ανάστροφου εκκρεμούς με παρουσία μεγάλου θορύβου δεν ισορροπούσε • Με μικρή μείωση του αρχικού θορύβου επιτεύχθηκε ο στόχος του να αντεπεξέρχεται το σύστημα κάτω από απρόσμενες συνθήκες • Το ανάστροφο εκκρεμές ισορροπεί μέχρι τις τιμές θορύβου που είναι 10-4, 10-5, 10-4 και 10-3 για τις εισόδους της «Θέσης του Βαγονιού», της «Γωνίας», της «Ταχύτητας του Βαγονιού» και της «Γωνιακής Ταχύτητας» αντίστοιχα

32. ΑΛΛΑΓΕΣ ΣΤΙΣ ΠΑΡΑΜΕΤΡΟΥΣ ΤΟΥ ΣΥΣΤΗΜΑΤΟΣ • Αλλάχθηκαν οι τιμές των παραμέτρων της μάζας του ανάστροφου εκκρεμούς (m) και του μήκους (lp)της ράβδου για φανεί πόσο αποδοτικός είναι ο ασαφής ελεγκτής. • Σκοπός ήταν να βρεθούν οι μέγιστες τιμές που μπορούσαν να δεχθούν και μέχρι τις οποίες θα ισορροπούσε το σύστημα του ανάστροφου εκκρεμούς. • Αρχικές τιμές τους ήταν m= 0,872 kg και lp= 0,5 m

33. ΑΛΛΑΓΗ ΣΤΗ ΜΑΖΑ ΤΟΥ ΑΝΑΣΤΡΟΦΟΥ ΕΚΚΡΕΜΟΥΣ • Αρχική τιμή ήταν 0,872 kg • Η ανώτερη τιμή που δόθηκε και το ανάστροφο εκκρεμές ισορροπούσε ήταν μέχρι 2,2 kg • Στα διαγράμματα διακρίνονται η θέση του βαγονιού (πάνω) και η γωνία του εκκρεμούς (κάτω) για τιμή της m=2,2 kg

34. ΑΛΛΑΓΗ ΣΤΟ ΜΗΚΟΣ ΤΗΣ ΡΑΒΔΟΥ ΤΟΥ ΑΝΑΣΤΡΟΦΟΥ ΕΚΚΡΕΜΟΥΣ • Αρχική τιμή της ήταν 0,5 m. • Η ανώτερη τιμή που δόθηκε και το ανάστροφο εκκρεμές ισορροπούσε ήταν μέχρι 1,2 m. • Στα διαγράμματα διακρίνονται η θέση του βαγονιού (πάνω) και η γωνία του εκκρεμούς (κάτω) για τιμή της lp=1,2 m

35. ΑΛΛΑΓΕΣ ΣΤΙΣ ΠΑΡΑΜΕΤΡΟΥΣ ΤΗΣ ΜΑΖΑΣ ΚΑΙ ΤΟΥ ΜΗΚΟΥΣ ΤΗΣ ΡΑΒΔΟΥ ΤΟΥ ΑΝΑΣΤΡΟΦΟΥ ΕΚΚΡΕΜΟΥΣ • Αρχικές τιμές τους ήταν m= 0,872 kg και lp=0,5 m • Οι μέγιστες τιμές που δόθηκαν στις δύο παραμέτρους και το ανάστροφο εκκρεμές ισορροπούσε ήταν m=1,8kgκαι lp=0,9 m. • Στα διαγράμματα διακρίνονται η θέση του βαγονιού (πάνω) και η γωνία του εκκρεμούς (κάτω) για τις τιμές αντίστοιχα m=1,8kgκαι lp=0,9 m.

36. ΑΠΟΤΕΛΕΣΜΑΤΑ ΔΟΚΙΜΩΝ ΑΠΟ ΤΙΣ ΑΛΛΑΓΕΣ ΤΩΝ ΠΑΡΑΜΕΤΡΩΝ

37. ΣΥΜΠΕΡΑΣΜΑΤΑ ΑΛΛΑΓΗΣ ΤΩΝ ΠΑΡΑΜΕΤΡΩΝ • Οι τιμές που μπορούν να πάρουν οι παράμετροι του συστήματος για να ισορροπεί το ανάστροφο εκκρεμές είναι: • Η μάζα ανάμεσα στις 0,872 – 2,2 kg • Το μήκος της ράβδου ανάμεσα στις0,5 - 1,2 m • Η μάζα του και το μήκος της ράβδου όταν αλλάζουν συγχρόνως, ανάμεσα στις 0,872 – 1,8 kg όσον αφορά τη μάζα του ανάστροφου εκκρεμούς και από το 0,5-0,9 m όσον αφορά το μήκος της ράβδου

38. ΑΝΑΚΕΦΑΛΑΙΩΣΗ • Σκοπός της εργασίας ήταν η εξισορρόπηση του ανάστροφου εκκρεμούς με την εφαρμογή ενός ευφυή ελεγκτή • Για την επίλυση του προβλήματος χρησιμοποιήθηκε το ασαφές μοντέλο Mamdani • Οι κανόνες του συστήματος δημιουργήθηκαν από την προσωπική εμπειρία του σχεδιαστή του συστήματος • Εισάγοντας θόρυβο στο σύστημα καθώς και αλλάζοντας τις παραμέτρους του, ελέγχθηκε η ευρωστία του ασαφή ελεγκτή

39. ΕΥΧΑΡΙΣΤΙΕΣ • Προς τον καθηγητή Δρ. Βολογιαννίδη για την πολύτιμη βοήθειά του • Προς τους καθηγητές του Τμήματος για τις γνώσεις που μας προσέφεραν κατά τη διάρκεια των σπουδών μας στο Τμήμα Πληροφορικής και Επικοινωνιών του ΤΕΙ Σερρών

40. ΕΥΧΑΡΙΣΤΩ ΠΟΛΥΓΙΑ ΤΗΝ ΠΡΟΣΟΧΗ ΣΑΣ