Metode de Optimizare de Tip Evolutiv Pentru Sistemele de Mari Dimensiuni

LUCRARE DE LICENȚĂ

Metode de optimizare de tip evolutiv pentru sistemele de mari dimensiuni

REZUMAT

1. INTRODUCERE

1.1. PROBLEMA DIMENSIUNII ÎN REZOLVAREA EFECTIVĂ A PROBLEMELOR DE OPTIMIZARE

1.2. CLASIFICAREA METODELOR DE REZOLVARE A PROGRAMELOR MARI

2. TEHNICI DE DESCOMPUNERE SAU DE PARTIȚIONARE

2.1. FOLOSIREA METODELOR DE DESCOMPUNERE

2.1.1. Scurtă prezentare

2.1.2. Exemplu de descompunere a unui sistem

2.1.3. Principiul de descompunere Dantzig-Wolfe

2.1.4. O interpretare economică a principiului de descompunere

2.2. METODE DE PARTIȚIONARE ȘI RELAXARE

2.2.1. Metode de partiționare și relaxare în programarea liniară

2.2.2. Proceduri de partiționare în programarea neliniară

2.3. PROGRAME NELINIARE CU VARIABILE DE CUPLARE (CAZ PARTICULAR)

2.3.1. Prezentare

2.3.2. Exemplu de abordare a unei probleme neliniare cu variabile de cuplare

3. METODE DE CĂUTARE EVOLUTIVE

3.1. METODE DE CĂUTARE DIRECTĂ

3.2. METODA NELDER—MEAD

3.2.1. Considerente teoretice

3.2.2. Aplicație la metoda Neder – Mead

3.3. METODA BOX

3.3.1. Considerente teoretic

3.3.2. Aplicație la Metoda Box (Complex)

3.4. ALGORITMULUI PSO (PARTICLE SWARM OPTIMISATION)

3.4.1. Considerente teoretice

3.4.2. Aplicație la metoda PSO

3.4.3. Concluzii pentru algoritmul PSO

3.4.4. Aplicație: Folosire ToolBox Matlab PSO

4. SISCON – PRODUS SOFTWARE PENTRU OPTIMIZARE

4.1. DESCRIERE SOFT (MANUAL DE UTILIZARE)

4.2. METODE DE CĂUTARE EVOLUTIVE

4.2.1. Aplicație la Metoda Nelder-Mead

4.2.2. Aplicație la Metoda BOX

4.3. TRANSFORMĂRI ECHIVALENTE PENTRU PROBLEME DE MARI DIMENSIUNI

4.3.1. Tehnici de descompunere

4.3.2. Tehnici de penalizare

4.3.3. Transformarea duală

5. STUDIU COMPARATIV AL ALGORITMILOR EVOLUTIVI

5.1. EXEMPLE

5.2. OBSERVAȚII

6. CONCLUZII

7. ANEXE – COD ALGORITMI

7.1. METODA NELDER – MEAD

7.2. METODA BOX

7.3. PSO (PARTICLE SWARM OPTIMISATION)

7.4. TOOLBOX MATLAB PSO

BIBLIOGRAFIE

Rezumat

În această lucrare este descrisă o anumită clasă de metode de optimizare și anume, metode de optimizare de tip evolutiv pentru sistemele de mari dimensiuni.

În Capitolul I este prezentată problema dimensiunii mari în rezolvarea efectivă a aplicațiilor de optimizare și sunt clasificate metodele de rezolvare ale acestor probleme.

Rezolvarea problemelor de mari dimensiuni se poate face în doua etape care vor fi tratate în capitolele următoare:

Aplicarea de tehnici de descompunere sau de partiționare în urma cărora rezultă probleme mai mici dar interconectate

Aplicarea de metode specifice problemelor cu restricții precum metodele de căutare evolutive

Capitolul al II-lea este dedicat tehnicilor de descompunere sau de partiționare. În acest capitol sunt prezentate:

modalitatea de folosire a metodelor de descompunere, se arată cum trebuie descompus un sistem pe un exemplu dat, este făcut un rezumat asupra principiului de descompunere Danzing-Wolfe și este dată o interpretare economică a acestui principiu

metode de partiționare și relaxare atât pentru programarea liniară cât și pentru cea neliniara, sunt explicate noțiunile de relaxare și restrângere iar în final este prezentat algoritmul lui Rosen pentru programe neliniare

programe neliniare cu variabile de cuplare (caz particular). Pentru a vedea cum o astfel de problemă conduce la un program neliniar; s-a considerat un exemplu simplificat cu două rafinării R1 și R2 care fabrică produsele P1 și P2, pentru a fi livrate la stațiile terminale T1 și T2 .

Capitolul al III-lea cuprinde metode de căutare directă. Focalizarea este făcută asupra metodelor de căutare directă multivariabile evolutive și sunt prezentate:

Metoda NELDER-MEAD

Metoda BOX

Algoritmul PSO (Particle Swarm Optimisation)

Pentru fiecare metodă în parte s-a realizat o scurtă prezentare din punct de vedere teoretic și au fost prezentate aplicații sugestive pentru a se arăta modul de evoluție și particularitățile algoritmilor

Capitolul al IV-lea cuprinde o prezentare produsului software pentru optimizare, SISCON. Pachetul de programe SISCON, este o aplicație destinată optimizării și conducerii eficiente a sistemelor, dezvoltată în C++ (Builder). Rolul meu a fost de a folosi acest program și de a arăta pe câteva exemple care sunt rezultatele obținute și care este modul lui de funcționare.

Capitolul al V-lea conține un studiu comparativ pentru cele trei metode prezentate anterior realizat pe mai multe exemple și sintetizate rezultate și informații legate de anumiți indicatori precum:

timp de rulare

memorie ocupată de algoritm

iterații până la atingerea optimului

eroare

importanța condițiilor inițiale

Capitolul al VI-lea este destinat concluziilor.

1. INTRODUCERE

1.1. Problema dimensiunii în rezolvarea efectivă a problemelor de optimizare

Principala cauză generatoare de dificultăți în rezolvarea problemelor de optimizare reale este dimensiunea: o anumită problemă este pur și simplu "prea mare". În programarea matematică, mărimea unei probleme este o chestiune relativă, depinzând de mulți parametri cum ar fi:

numărul de variabile și numărul de restricții;

complexitatea expresiilor funcției obiectiv și a restricțiilor;

performanțele echipamentului de calcul: hardware și software.

Ceea ce înseamnă „mare" depinde de capacitatea echipamentului de calcul disponibil. Totuși, odată cu folosirea pe scară mai intensivă a programării matematice, au fost formulate probleme care suprasolicită chiar cei mai buni algoritmi și calculatoarele cele mai mari. Totuși, au fost rezolvate probleme având 104 ecuații și 106 variabile. Pentru realizarea unor astfel de performanțe este necesară construirea unor algoritmi speciali care să utilizeze structura sistemului.

Utilizarea modelelor matematice în studiul unor situații reale a condus la programe matematice care sunt greu de rezolvat și pe cele mai puternice calculatoare. Din fericire, marea majoritate a problemelor de optimizare "mari" au o "structură specială" care, în programarea liniară de exemplu, înseamnă:

densitate mică (câteva procente) a constantelor numerice nenule din matricea coeficienților;

gruparea elementelor nenule în blocuri așezate "pe diagonală", exceptând un număr mic de linii și coloane ;

număr foarte mare de variabile și relativ puține restricții sau invers, multe restricții și puține variabile ceea ce este echivalent cu a avea multe linii sau multe coloane, dar nu amândouă în același timp

mulțimea coloanelor (liniilor) este definită prin anumite relații explicite, de exemplu inegalități liniare.

Această situație permite căutarea eficientă a unei linii sau coloane oarecare, elementul dorit fiind generat la cerere.

Trebuie spus că dacă un program liniar mare nu are o structură specială, sarcina colectării datelor este practic aproape imposibil de realizat; astfel pentru un program cu 104 restricții și 106 variabile, matricea coeficienților ar avea 1010 intrări și în cazul unei densități de 100% ar necesita 1010 date numerice de adunat, sortat și prelucrat!!

În cazul neliniar programele se clasifică mai greu. Se poate defini o matrice al cărei element (i, j) este egal cu unu dacă variabila i apare în restricția j și zero în rest, și această matrice va reflecta structura problemei. Are de asemenea importanță să se știe care variabile apar liniar și dacă unele funcții neliniare sunt separabile. Apar în mod obișnuit problemele hibride, cum ar fi de exemplu problemele care devin liniare cu structură specială atunci când un număr relativ mic de variabile neliniare se fixează.

1.2. Clasificarea metodelor de rezolvare a programelor mari

Metodele de rezolvare a programelor matematice mari se pot împărți în două clase:

tehnici de descompunere sau de partiționare

metode directe

Metodele directe specializează o procedură generală adaptând-o la specificul unei anumite clase de probleme de optimizare.

În categoria metodelor de căutare directă multivariabile amintim:

metode de căutare multivariabile liniare

metode de căutare evolutive

metode de gradient

Exemplul tipic îl constituie algoritmul simplex; este știut că principala problemă de calcul care apare în aplicarea lui o constituie modul de manipulare a inversei bazei curente. În cazul unei "structuri speciale" este posibil ca dimensiunea acestei matrice să se reducă semnificativ. . Printre metodele care procedează în acest fel amintim procedurile pentru probleme cu variabile mărginite superior și metodele de triangularizare compactă a bazei.

Să considerăm cazul unui program liniar cu variabile superior mărginite:

Abordarea "clasică" presupunea transformarea condițiilor de limitare superioară în egalități:

Rezulta un program cu m+n restricții și 2n variabile, ale cărui baze erau matrice de ordinul m+n.

Totuși forma particulară a restricțiilor de limitare superioară a putut fi exploatată eficient într-o specializare a algoritmului simplex în care inversa bazei curente are dimensiunea egală cu numărul m al restricțiilor propriu zise.

Metodele indirecte, bazate pe descompunerea problemei mari în subprobleme mai mici, interconectate. Subproblemele pot fi rezolvate independent (și dacă este posibil chiar simultan) dar faptul că ele interacționează implică existența unui mecanism (problemă) de coordonare. Astfel, rezolvarea problemei originale mari se face "la două niveluri":

la primul nivel – inferior – se rezolvă subproblemele; rezultatele sunt "comunicate"

la al doilea nivel – superior – care "analizează" aceste rezultate și transmite nivelului inferior noi parametri.

La nivelul unu are loc o reluare a calculelor (reoptimizare); noile rezultate sunt trimise nivelului superior care le analizează ș.a.m.d.

Important este faptul că acest proces iterativ este convergent în sensul că într-un număr finit de pași (de dialoguri între cele două niveluri), nivelul coordonator "anunță" găsirea soluției optime.

Deoarece subsistemele interacționează, rezolvarea subproblemelor nu va furniza în general soluția corectă. Abordarea la mai multe niveluri a problemei propune ca ținând seama de interacțiuni să se definească unul sau mai multe subsisteme ca fiind situate la „nivelul doi", care să influențeze într-un anumit mod subsistemele originale, definite ca fiind la nivelul întâi. Această influență poate lua forme numeroase în funcție de problema originală, de natura descompunerii de la primul nivel și de aceasta trebuie neapărat să se țină cont atunci când inițial se descompune sistemul. Scopul nivelului al doilea este de a coordona în așa fel acțiunile unităților de la primul nivel, încât să se obțină soluția problemei originale. Toate structurile organizatorice mari operează în acest mod.

Această idee poate fi extinsă prin definirea subsistemelor de la nivelul trei, fiecare dintre acestea coordonând un număr de unități de la nivelul doi etc. Structura rezultată are forma unei piramide a unităților decizionale, cum este arătat în Figura 1 – 1.

Figura 1 – 1 Structura de control la mai multe niveluri

Toți algoritmii de descompunere existenți au structura la două niveluri. Coordonatorul influențează subproblemele variind parametrii din funcțiile lor obiectiv, adăugind restricții etc. Avantajele acestei abordări includ libertatea de a rezolva fiecare subprogram alegând la dorință orice algoritm, precum și posibilitatea de a scrie și corecta programe de calculator independente mai mici și de a reduce cantitatea de memorie cerută. Principalul dezavantaj constă în faptul că aceste probleme trebuie rezolvate de un anumit număr de ori, astfel că trebuie utilizat mult timp calculator.

Astfel, în paragrafele următoare, vom considera un număr de exemple. Deseori va fi util să se facă abstracție de semnificația fizică a acestora și să le discutăm folosind aceeași terminologie. Conceptele de analiză a activității care furnizează o bază sistematică pentru construirea de modele de programare matematică permit acest lucru.

2. Tehnici de descompunere sau de partiționare

2.1. Folosirea metodelor de descompunere

2.1l doi", care să influențeze într-un anumit mod subsistemele originale, definite ca fiind la nivelul întâi. Această influență poate lua forme numeroase în funcție de problema originală, de natura descompunerii de la primul nivel și de aceasta trebuie neapărat să se țină cont atunci când inițial se descompune sistemul. Scopul nivelului al doilea este de a coordona în așa fel acțiunile unităților de la primul nivel, încât să se obțină soluția problemei originale. Toate structurile organizatorice mari operează în acest mod.

Această idee poate fi extinsă prin definirea subsistemelor de la nivelul trei, fiecare dintre acestea coordonând un număr de unități de la nivelul doi etc. Structura rezultată are forma unei piramide a unităților decizionale, cum este arătat în Figura 1 – 1.

Figura 1 – 1 Structura de control la mai multe niveluri

Toți algoritmii de descompunere existenți au structura la două niveluri. Coordonatorul influențează subproblemele variind parametrii din funcțiile lor obiectiv, adăugind restricții etc. Avantajele acestei abordări includ libertatea de a rezolva fiecare subprogram alegând la dorință orice algoritm, precum și posibilitatea de a scrie și corecta programe de calculator independente mai mici și de a reduce cantitatea de memorie cerută. Principalul dezavantaj constă în faptul că aceste probleme trebuie rezolvate de un anumit număr de ori, astfel că trebuie utilizat mult timp calculator.

Astfel, în paragrafele următoare, vom considera un număr de exemple. Deseori va fi util să se facă abstracție de semnificația fizică a acestora și să le discutăm folosind aceeași terminologie. Conceptele de analiză a activității care furnizează o bază sistematică pentru construirea de modele de programare matematică permit acest lucru.

2. Tehnici de descompunere sau de partiționare

2.1. Folosirea metodelor de descompunere

2.1.1. Scurtă prezentare

În ultimul deceniu a căpătat o extindere din ce în ce mai largă o tehnică de optimizare bazată pe descompunerea unor probleme de optimizare pentru instalații sau sisteme cu structuri mari și de complexitate ridicată în probleme optimizare ale sistemelor mai simple ale instalațiilor respective; această tehnică este denumită „tehnică de descompunere" sau „tehnică de decompoziție".

O asemenea descompunere poate prezenta avantaje importante pentru obținerea soluțiilor optimizării, întrucât problemele de optim ale subsistemelor se rezolvă mai ușor, datorită numărului redus de variabile care intervin în fiecare subsistem; pentru obținerea unui optim pe ansamblu, pe lângă optimizarea funcționării fiecărui subsistem este necesară și o acțiune de coordonare centrală, pentru compensarea legăturilor de interacțiune dintre subsisteme.

Se obține astfel o structură ierarhizată, cu două sau mai multe niveluri; în cazul a două niveluri, menționat anterior, la nivelul inferior are loc optimizarea funcționării fiecărui subsistem, iar la nivelul superior se efectuează acțiunea de coordonare pentru obținerea optimului pe ansamblu. Considerentele referitoare la ierarhizarea pe două niveluri se pot generaliza ușor pentru cazul mai multor niveluri

.

Figura 2 – 1

Structura unui sistem cu. două niveluri este reprezentată în Figura 2 – 1. La nivelul inferior se găsesc subsistemele 1, 2, 3 care asigură optime locale, de exemplu prin maximizarea unor criterii , iar nivelul superior este constituit de coordonatorul 4, care asigură un optim global (schimbând informații în ambele sensuri cu subsistemele 1, 2, 3, informații reprezentate simbolic prin săgeți) de exemplu prin minimizarea unui criteriu .

2.1.2. Exemplu de descompunere a unui sistem

Pentru a ilustra descompunerea unui sistem complex în două subsisteme, în Figura 2 – 2 este reprezentată o variantă de descompunere a unei instalații tehnologice . Cele două subsisteme formate, subsistemul 1 și subsistemul 2, sunt legate prin două mărimi: debitul de intrare în coloana C, notat cu la ieșirea din subsistemul 1 și cu la intrarea în subsistemul 2, și fracțiunea recirculată din produsul de bază, notată cu , la ieșirea din subsistemul 2 și cu , la intrarea în subsistemul 1. Blocul de coordonare, de la nivelul superior, nu este reprezentat în Figura 2 – 2.

La restricțiile anterioare de tip egalitate se adaugă alte restricții. Una dintre ele rezultă din ecuația bilanțului de material pentru debitul care iese din subsistemul 1:

(2-1)

această ecuație conducând la restricția:

(2-2)

Figura 2 – 2

Alte două restricții de tip egalitate rezultă din necesitatea egalității debitelor de ieșire din unul dintre subsisteme și de intrare în celălalt, obținându-se egalitățile:

și (2-3)

respectiv restricțiile

și (2-4)

Pentru subsistemul 2 rezultă o singură ecuație de bilanț de material:

(2-5)

din aceasta obținându-se restricția

(2-6)

Alegerea unor proceduri de coordonare eficiente, la nivelul superior, are o importanță deosebită pentru reducerea efortului de calcul pentru optimizarea prin metode de descompunere.

Conceptul de „structură multinivel" nu este singurul care intervine în cadrul optimizărilor ierarhizate. Astfel, pe lângă structura multinivel (în limba engleză, multilevel structure) este utilizată și „structura multistrat" sau „multialgoritm" (în limba engleză, multilayer structure) pentru realizarea optimizărilor ierarhizate.

Structura multistrat realizează, de asemenea, o ierarhizare și anume pe „straturi" sau „algoritmi", la diversele straturi intervenind algoritmi cu scopuri diferite și care acționează la intervale de timp diferite. Pentru ilustrare, în Figura 2 – 3 este prezentată o structură multistrat pentru automatizarea instalației tehnologice din blocul F, care include și elemente de execuție și traductoare.

Stratul 1 conține algoritmii de reglare automată și acționează în permanență pentru ca valorile mărimilor reglate, constituind vectorul , să urmărească valorile prescrise prin mărimile de referință, constituind vectorul ; în acest scop, stratul 1 transmite elementelor de execuție din blocul vectorul mărimilor de comandă .

Asupra blocului acționează și vectorul perturbărilor .

Figura 2 – 3

Stratul 2 include un algoritm de optimizare, care pe baza informațiilor și – primite de la blocul – și a informațiilor – primite de la mediul înconjurător – elaborează valorile optime prescrise pentru mărimile reglate , transmițând stratului 1 aceste valori sub forma vectorului mărimilor de referință ; valorile optime menționate sunt calculate de algoritmul de optimizare în ipoteza că parametrii instalației tehnologice (factori de amplificare, constante de timp, timp mort, constituind un vector ) nu suferă variații în timp.

Dacă această ipoteză nu este valabilă și parametrii instalației tehnologice suferă variații sub influența unor perturbări parametrice, atunci devine necesară prezența unui al treilea strat, respectiv stratul 3 din Figura 2 – 3, care realizează un algoritm de adaptare. Pe baza informațiilor primite de la blocul și a informațiilor primite de la mediul înconjurător, algoritmul 3 determină valorile parametrilor și le transmite stratului 2, pentru ca algoritmul de optimizare din acest strat să ia în considerare valorile respective, uneori fiind modificată în mod corespunzător structura acestui algoritm.

Stratul 3 asigură astfel o actualizare a valorilor parametrilor modelului matematic al instalației tehnologice.

În unele cazuri, ca urmare a primirii valorilor stratul 2 poate transmite stratului 1 și alte semnale, pe lângă vectorul .

Pe lângă cele trei straturi menționate poate exista și un al patrulea strat, cu algoritm de auto-organizare, care în funcție de anumiți factori exteriori poate asigura eventuale modificări ale structurii sistemului ierarhizat.

În sistemele cu structură multistrat fiecare strat intervine la intervale de timp diferite, intervalele fiind-cu atât mai mari, cu cât stratul are o poziție mai înaltă. De asemenea, fiecare strat folosește un grad de simplificare diferit pentru modelul instalației tehnologice, modelul fiind cu atât mai simplificat și mai abstractizat cu cât poziția stratului este mai înaltă.

În unele cazuri structura multistrat are un aspect preponderent temporal. Astfel, de exemplu în cazul optimizării utilizării resurselor de apă, stratul superior determină soluția optimă pe termen lung, de un an, folosind un model simplificat, în care toate rezervoarele mai puțin importante sunt înlocuite printr-un rezervor echivalent. Pe baza calculelor de optimizare din stratul superior se transmite stratului mijlociu o soluție care trebuie realizată la sfârșitul primei luni a anului, soluție care nu este însă detailată, întrucât rezervoarele de dimensiuni medii și mici nu au fost considerate la volumele lor reale

.

Stratul mijlociu calculează mai precis soluția de optimizare pe termen de o lună, considerând datele exacte ale rezervoarelor medii și înlocuind toate rezervoarele mici printr-un rezervor echivalent, rezultatul optimizării pentru sfârșitul primei zi a lunii fiind transmis stratului inferior. în cadrul stratului inferior se calculează soluția detailată pe termen de o zi – ținând seama de dimensiunile reale ale rezervoarelor mici – și se asigură realizarea acestei soluții prin intermediul echipamentelor de reglare automată ale instalațiilor de distribuție a apei.

Pe măsura realizării soluțiilor calculate (deci după fiecare zi, respectiv după fiecare lună) se transmit semnale de reacție de la stratul inferior la cel mediu și de la cel mediu la cel superior, aceste informații fiind utilizate ca valori inițiale – de straturile care le primesc – pentru recalcularea soluțiilor optime pe intervalele aferente, de o lună și respectiv un an.

2.1.3. Principiul de descompunere Dantzig-Wolfe

2.1.3.1. Prezentare

Publicarea în 1960 de către Dantzig și Wolfe a principiului de descompunere a inițiat în perioada care a urmat o activitate intensivă în domeniul programării matematice de dimensiuni mari. Procedura este în special eficientă atunci când se aplică la programe liniare ale căror matrice ale coeficienților au structură unghiulară, adică formează unul sau mai multe blocuri independente legate prin ecuații de cuplare. Ea operează prin formarea unui „program principal" echivalent, care are un număr de linii numai cu puțin mai mare decât numărul ecuațiilor de cuplare din problema generală, dar care are un număr foarte mare de coloane.

Acest program se rezolvă fără a tabela toate aceste coloane, generându-le însă ori de câte ori metoda simplex le necesită, tehnică pe care o vom numi „generare de coloane". Algoritmul care rezultă implică o iterație între o mulțime de subprobleme independente ale căror funcții obiectiv conțin parametri variabili și programul principal. Sub-problemele primesc o mulțime de parametri de la programul principal. Ele își trimit soluțiile către programul principal, care le combină cu soluțiile anterioare într-un mod optimal și calculează noi prețuri. Acestea sunt din nou trimise subproblemelor și iterația continuă până când este verificat un test de optimalitate.

Procedura are o interpretare economică elegantă, în care programul principal coordonează acțiunile subproblemelor punând prețuri pe resursele folosite de ele.

2.1.3.2. Expunerea principiului de descompunere

Se consideră un program liniar a cărui matrice a coeficienților are o structură unghiulară în – blocuri, cu :

(2-8)

Orice program liniar poate fi considerat că are această formă cu , obținută prin partiționarea restricțiilor în două submulțimi :

să se minimizeze (2-9)

în raport cu

(2-10)

(2-11)

(2-12)

Se presupune că poliedronul convex

(2-13)

este mărginit; această ipoteză va fi ulterior slăbită. Atunci, orice element al lui este scris

, (2-14)

unde

, (2-15)

iar sunt puncte extreme ale poliedronului .

Se obține

(2-16)

Substituind (2-14) în (2-9), se obține exprimarea lui în termenii variabilelor :

(2-17)

Definind

(2-18)

, (2-19)

relațiile (2-15) – (2-17) conduc la un program liniar în :

să se minimizeze , (2-20)

în raport cu

, (2-21)

, (2-22)

, (2-23)

Acest program, numit program principal, este echivalent cu (2-18) cel original. El are numai linii. Față de problema originală câre are linii, acesta reprezintă o restrângere a dimensiunii în cazul în care este mare. De asemenea, el are atâtea coloane câte puncte extreme are poliedronul , care pot fi în număr de mii dacă este mare.

În loc de a tabela toate aceste coloane folosim o tehnică de generare de coloane, producând coloane care să intre în bază ori de câte ori este nevoie de ele. Pentru a vedea cum se face acest lucru se consideră coeficientul de cost relativ pentru variabila :

. (2-24)

Se partiționează astfel:

, (2-25)

unde corespunde restricțiilor (2-21) iar scalarul restricției (2-22). Atunci, folosind definițiile lui și p din (2-18) și (2-19), poate fi scris în modul următor :

(2-26)

Conform criteriului simplex, pentru a alege variabila care intră în bază trebuie să găsim:

(2-27)

Reamintim că este un punct extrem al lui și observăm că este liniar în . Deoarece o soluție optimă a unui program liniar a cărui mulțime restricție este mărginită este întotdeauna un punct extrem al acestei mulțimi, minimizarea din (2-27) este echivalentă cu subproblema

Să se minimizeze (2-28)

în raport cu

(2-29)

Pentru a găsi o coloană care să intre în baza programului principal se rezolvă această subproblemă obținând o soluție . Atunci coloana care intră în bază este

(2-30)

coeficient de cost fiind (2-31)

Această procedură devine deosebit de atractivă dacă , numărul de blocuri independente din structura unghiulară, este mai mare decât unu, adică dacă problema care trebuie rezolvată are forma

să se minimizeze (2-32)

în raport cu

,

,

, (2-33)

:

(2-34)

Atunci subproblema (2-28) – (2-29) devine

să se minimizeze , (2-35)

în raport cu (2-36)

Deoarece obiectivul (2-35) este aditiv separabil în și restricțiile (2-36) asupra lui sunt independente, această problemă se reduce la subprobleme independente :

să se minimizeze (2-36)

în raport cu (2-37)

Algoritmul de descompunere. Vom putea acum formula un algoritm la două niveluri pentru soluționarea problemei (2-32) – (2-34), cu programul principal la nivelul secund și cu subproblemele (2-36) – (2-37) la primul nivel. Se presupune că avem la dispoziție o soluție admisibilă de bază inițială pentru programul principal (2-20) – (2-23) cu matricea bazică și cu multiplicatorii simplex .

Pasul 1. Folosind multiplicatorii simplex se rezolvă subproblemele(2-36) – (2-37) obținând soluțiile și valorile optimale ale obiectivului .

Fie '. (2-38)

Pasul 2. Se calculează min care, conform cu (2-26), este

(2-39)

Dacă stop (2-40)

Soluția optimala pentru (2-32) – (2-34) este

, (2-41)

unde sunt punctele extreme ale lui corespunzătoare valorilor bazice.

Pasul 3. Dacă min < 0, se formează coloana

(2-42)

Aceasta se transformă prin multiplicare cu și folosind operația uzuală de pivotare simplex se obține o nouă inversă de bază și un nou vector al multiplicatorilor simplex. Se revine la pasul 1 și se repetă.

Dacă programul principal este nedegenerat (sau – perturbat), atunci fiecare iterație îl va descrește pe cu o cantitate diferită de zero. Deoarece există doar un număr finit de baze posibile și nici una nu se repetă, principiul de descompunere va găsi soluția optimală într-un număr finit de iterații.

Remarcăm că soluția optimală nu este neapărat una dintre propunerile subproblemelor. Din (2-41) rezultă că este o anumită combinație convexă a unui număr de astfel de soluții. Astfel, programul principal nu poate obține optimul global pur și simplu prin trimiterea, prețurilor corespunzătoare către subprobleme; el trebuie să aibă libertate să combine soluțiile subproblemelor într-un plan general. În acest sens descompunerea obținută aici nu poate fi considerată ca o descentralizare completă a procesului de luare a deciziei. O denumire mai potrivită dată de Dantzig acestei descompuneri este „planificare centralizată fără informație completă la centru".

2.1.4. O interpretare economică a principiului de descompunere

Să considerăm o economie cu mai mulți agenți. Fiecare agent operează un număr de activități de pe urma cărora obține un venit. Operarea activităților implică utilizarea unor resurse disponibile în cantități limitate. Firește, obiectivul fiecărui agent este maximizarea venitului propriu.

Este clar că dacă fiecare agent ar deține controlul asupra tuturor resurselor necesare lui atunci maximizarea venitului la scara întregii economii s-ar reduce la maximizarea venitului fiecărui agent în parte.

În realitate lucrurile stau altfel.

Fiecare agent deține controlul asupra anumitor resurse: capacități proprii de producție, forța de muncă angajată, resurse financiare proprii, unele materii prime utilizate în exclusivitate. Acestea vor fi numite resurse specifice.

Pe lângă resursele specifice, fiecare agent utilizează și alte resurse care nu sunt la dispoziția sa exclusivă; aceste resurse sunt procurate de pe piață, la concurență cu ceilalți agenți, datorită faptului că sunt disponibile în cantități limitate. Acestea vor fi denumite resurse comune.

În acest context, problema principială care se pune este de a stabili cum vor fi repartizate resursele comune între agenți astfel încât, la scara întregii economii, venitul să fie maxim.

Într-o economie centralizată, repartizarea resurselor comune este făcută de stat care indică fiecărui agent ce și cât să producă.

Ne propunem să arătăm cum se face repartizarea comune într-o economie descentralizată în care autoritatea centrală nu mai deține controlul asupra acțiunilor agenților.

Ne vom situa în cazul ideal al unei economii liniare, caracterizate prin următoarele ipoteze:

pentru fiecare activitate, consumurile de resurse și venitul sunt direct proporționale cu nivelul la care este operată activitatea;

nivelul de operare al unei activități poate fi reprezentat de orice număr real nenegativ;

veniturile agenților nu se condiționează reciproc și sunt egale cu suma veniturilor activităților fiecăruia. Venitul la scara întregii economii este suma veniturilor agenților.

Pentru simplitate vom presupune că în economia studiată există numai doi agenți.

Introducem notațiile:

= vectorii (coloană) nivelelor de operare ale activităților celor doi agenți;

= vectorii (coloană) cantităților disponibile din resursele specifice;

= matricele consumurilor unitare de resurse specifice;

= metricile consumurilor unitare din resursele comune;

= vectorul (coloană) cantităților disponibile de resurse comune;

= vectorii (linie) veniturilor unitare corespunzătoare celor doi agenți.

Evident, nivelele de operare ale activităților agenților sunt condiționate de disponibilele de resurse specifice:

(2-43)

și în plus:

(2-44)

Vectorii care satisfac se vor numi programe posibile de activitate.

Un cuplu de programe posibile () devine realizabil numai dacă necesarul de resurse comune se încadrează în disponibilul dat adică:

(2-45)

Venitul total pe economie are expresia:

(2-46)

Reunind expresiile de mai sus se obține următorul program liniar:

(2-47)

care modelează problema repartizării resurselor comune în vederea maximizării venitului pe întreaga economie.

Observăm că matricea programului are structura bloc diagonală cu restricții de cuplare, identică cu structura pe care se poate aplica principiul de descompunere Dantzig – Wolfe.

Din punct de vedere formal, problema repartizării resurselor comune într-o economie descentralizată înseamnă rezolvarea programului (P) în condițiile în care nici agenții nici autoritatea centrală nu au informații complete asupra acestuia. Astfel:

agentul 1 controlează (cunoaște) ;

agentul 2 controlează (cunoaște) ;

autoritatea centrală controlează (cunoaște) .

Maximizarea venitului fiecărui agent, ținând seama numai de resursele sale specifice, revine formal la rezolvarea programelor:

(2-48)

dar nu rezolvă problema repartizării resurselor comune deoarece, dacă și sunt soluțiile optime ale programelor și , este posibil ca:

(2-49)

În continuare vom arăta – în principiu – cum poate fi rezolvat programul (P) în situația în care nici autoritatea centrală și nici agenții nu dețin informații complete asupra programului!

Vom presupune că:

între autoritatea centrală și agenți există o cooperare în sensul unui schimb de informații privind "intențiile"de acțiune;

autoritatea centrală își asumă rolul de arbitru în următorul sens: ea "anunță" un sistem de prețuri pe resursele comune iar agenții iau aceste prețuri ca date și își diminuează veniturile cu valoarea resurselor comune solicitate. Fie u vectorul (linie) al acestor prețuri. Atunci:

agentul 1, pentru susținerea unui program posibil (posibil ),va trebui să "plătească" valoarea astfel că venitul său "net" va fi:

(2-50)

analog, venitul agentului 2, rezultat din programul posibil() va fi:

(2-51)

Maximizarea acestor venituri nete înseamnă rezolvarea programelor modificate:

(2-52)

Agenții comunică autorității centrale propunerile optimale și . În principiu, autoritatea centrala analizează oportunitatea luării în considerare a acestor propuneri pentru maximizarea venitului pe economie și poate decide modificarea sistemului de prețuri, mărind prețurile resurselor comune "intens" solicitate.

Noile prețuri sunt comunicate agenților; aceștia vor căuta noi soluții care să le maximizeze veniturile nete "corectate". Evident, prin creșterea prețurilor la anumite resurse comune, cererile "excesive" din aceste resurse vor fi "temperate". Formal, cele spuse înseamnă reluarea programelor cu u modificat!

Important este că într-un număr finit de asemenea "dialoguri = schimburi de informații" între agenți și autoritatea centrală, vor rezulta soluțiile și care maximizează venitul total pe economie. În general, și nu coincid cu una sau alta dintre propunerile agenților ( propuneri făcute în cadrul dialogului sus amintit) ci sunt combinații convexe ale acestora.

Totodată va rezulta și un sistem de prețuri pe resursele comune în raport cu care și maximizează veniturile nete individuale ale agenților! Spunem că tripletul (,, ) reprezintă un echilibru pentru economia (liniară) considerată.

Comparând această schemă cu cea din Figura 2 – 4 se constată că cele spuse mai sus constituie o interpretare economică a principiului de descompunere Dantzig – Wolfe.

Dialogul dintre autoritatea centrală și agenți poate fi reprezentat astfel:

Figura 2 – 4

2.2. Metode de partiționare și relaxare

2.2.1. Metode de partiționare și relaxare în programarea liniară

2.2.1.1. Introducere

În continuare vor fi studiate programe liniare ale căror matrice au structură bloc diagonală, legătura dintre blocuri realizându-se prin variabile de cuplare sau prin variabile și restricții de cuplare. Fie un program cuplat în ambele moduri:

să se maximizeze (2-62)

în raport cu

Pentru rezolvarea acestei probleme se poate aplica principiul de descompunere Dantzig-Wolfe. Dacă problema are numai variabile de cuplare, atunci procedura Dantzig-Wolfe poate fi aplicată dualei (care are numai restricții de cuplare). Dacă sunt prezente atât variabile cit și restricții de cuplare, atunci, după cum se vede, matricea restricțiilor este partiționată în (2-62) și (2-63). Relațiile (2-64) devin restricții ale subproblemelor, iar (2-63) devine program principal.

Subproblema poate fi tratată cu ajutorul metodelor din acest capitol sau prin dualizare și apoi din nou prin aplicarea metodei Dantzig-Wolfe în oricare dintre cele două cazuri, rezultatul general este un algoritm, ,,la trei niveluri", care implică o iterație (care rezolvă subproblemă) situată în interiorul unei bucle iterative externe (care rezolvă programul principal). Deoarece iterația internă este destul de complexă, pot apărea dificultăți de calcul cauzate în special de faptul că metoda Dantzig-Wolfe converge deseori lent.

În continuare vom furniza mijloace utile de rezolvare a problemei (2-62) – (2-64). Ele se numesc metode de partiționare deoarece variabilele problemei se partiționează la început în două submulțimi, variabilele de cuplare și variabilele bloc . Variabilele sunt din nou partiționate în mulțimi dependente și independente prin specificarea matricelor de bază în interiorul fiecărei matrice . Această partiționare permite ca restricțiile bloc și variabilele dependente să fie eliminate din program, obținându-se o problem redusă, de dimensiuni mai mici.

Se folosește adjectivul relaxare, deoarece prin eliminarea variabilelor dependente se pierde controlul asupra restricțiilor de nenegativitate. Conceptul de relaxare a restricțiilor poate fi plasat într-un cadru unitar care permite o mai bună înțelegere a aplicării lui la aceste metode. Din această cauză vom discuta, pentru început, acest concept.

2.2.1.2. Relaxare

Se consideră programul concav general

să se maximizeze (2-65)

în raport cu problema (P) (2-66)

(2-67)

unde și sunt funcții concave de vectorul – dimensional , iar este o submulțime din .

Considerăm aici cazul în care restricțiile sunt suficient de numeroase pentru a produce dificultăți în rezolvare. Atunci o strategie rezonabilă de rezolvare constă în relaxarea (adică eliminarea temporară) unora dintre acest restricții și din rezolvarea programului ce le conține pe cele rămase.

Dacă această problemă este neadmisibilă, la fel va fi și cea originală. Dacă aceasta este admisibilă și dacă soluția rezultată satisface restricțiile relaxate, atunci soluția este optimală pentru problema P. În cazul în care restricțiile relaxate nu sunt satisfăcute impunem una sau mai multe dintre aceste restricții și repetăm procedura. Relaxarea apare ca o strategie efectivă atunci când se știe că la optim un număr relativ mic de restricții sunt active.

O astfel de situație este adevărată în cazul programelor liniare nedegenerate (exact restricții sunt active) și, în general, pentru programele neliniare manifestarea ei poate fi chiar mai pregnantă.

Geoffrion a formalizat această idee după cum urmează. Fie M mulțimea întregilor {1,..,m}; se notează cu o submulțime a lui și fie programul

să se maximizeze (2-68)

în raport cu , (2-69)

(2-70)

Presupunem că poate fi găsită o mulțime inițială astfel încât are un suprem finit și că toate supremurile finite sunt atinse. Atunci o strategie de relaxare care permite atât eliminarea cât și adăugarea de restricții decurge după cum urmează :

Pasul 0. Se ia și se alege o mulțime inițială . astfel încât este mărginită superior în raport cu restricțiile problemei .

Pasul 1. Se rezolvă . Dacă este neadmisibilă, atunci și este neadmisibilă. În caz contrar se obține o soluție optimală .

Pasul 2. Dacă , pentru , atunci soluția este optimală pentru problema originală .

Pasul 3. În caz contrar, fie o submulțime a lui ce conține cel puțin un indice al unei restricții nesatisfăcute și fie

, (2-71)

Pasul 4. Dacă , atunci se înlocuiește prin și se revine la Pasul 1.

Pasul 5. Dacă , atunci se înlocuiește prin , se înlocuiește prin și se revine la Pasul 1.

Această procedură adaugă una sau mai multe restricții încălcate și, dacă descrește, atunci elimină acele restricții din care nu sunt active în . Faptul că nu elimină nici o restricție în cazul în care garantează, după cum se va vedea, că numai un număr finit de probleme trebuie rezolvate.

Aspecte privind calculul. Pentru ca relaxarea să fie eficientă trebuie să avem la dispoziție proceduri efective pentru identificară restricțiilor nesatisfăcute (adică pentru alegerea mulțimii ) și pentru soluționarea problemelor succesive . Există o multitudine de alegeri ale mulțimii . Un mod obișnuit de alegere este acela al indicelui restricției celei mai nesatisfăcute, la fel ca în metoda simplex dual. Această alegere este folosită frecvent atunci când numărul de restricții este cu adevărat mare și când mulțimea funcțiilor restricție are o anumită formă specială.

Deseori, restricția cea mai nesatisfăcută poate fi generată fără evaluarea tuturor funcțiilor restricție, prin rezolvarea unei subprobleme convenabile. Această „generare de linii" este duală cu procedura generării da coloane discutată anterior. Astfel experiența de calcul privind diferitele tactici de evaluare a coloanelor în programele liniare se poala utiliza și aici.

Alegerea la întâmplare a unei singure restricții nesatisfăcute este probabil mult mai puțin efectivă decât alegerea celei mai nesatisfăcute restricții. De asemenea, dacă se alege astfel încât să conțină mai mulți indici, cum ar fi de exemplu indicii celor mai nesatisfăcute restricții, atunci modificarea valorii de la iterație la iterație va crește. Desigur, va crește și munca necesitată de soluționarea problemelor succesive .

O opțiune care este utilă atunci când numărul restricțiilor este mare constă în identificarea celor mai nesatisfăcute restricții, după care se aplica relaxarea unei probleme restrânse constând din aceste restricții plus cele curente din . În interiorul acestei probleme restrânse poate fi de exemplu ales ca fiind indicele celei mai nesatisfăcute restricții. În forma ei primală această tactică constă în găsirea a coloane cu cele mai mici costuri reduse negative și în efectuarea iterațiilor simplex standard asupra unei probleme în care acestea reprezintă coloanele nebazice inițiale.

Această tactică se numește „evaluare multiplă" sau „suboptimizare" și conduce la o reducere substanțială a timpului calculator atunci când se folosește forma produs a inversei.

Opțiunea de a impune mai mult decât o singură restricție nesatisfăcută se implementează ușor în cazul în care relaxarea se aplică problemelor liniare cu restricții de cuplare sau/și cu variabile de cuplare.

2.2.1.3. Restrângere

Este important de remarcat faptul că tactica de relaxare a unei restricții primele corespunde unei tactici de a restrânge o restricție duală ca să devină . Tactica de restrângere este utilă atunci când se aplică în primală, și poate fi formalizată în aproape același mod ca relaxarea. Pentru aceasta să definim subproblemă în felul următor :

să se maximizeze (2-72)

în raport cu

, (2-73)

unde din nou este o submulțime a mulțimii .

Pentru a asigura că oricare este concavă, trebuie să presupunem că toate funcțiile sunt liniare. Orice restricție neliniară se încorporează în mulțimea . Să presupunem că poate fi găsită o mulțime nevidă inițială așa încât este admisibilă, că toate supremurile finite sunt atinse și dacă are o soluție optimală, atunci are și un vector optimal de multiplicatori Lagrange. Atunci, procedura de restrângere decurge în felul următor :

Pasul 0. Se ia și se alege o mulțime inițială astfel încât este admisibilă.

Pasul 1. Se rezolvă . Dacă funcția obiectiv este nemărginită, problema originală este nemărginită. În caz contrar se obține o soluție optimală .

Pasul 2. Fie multiplicatorul Lagrange optimal asociat restricției . Dacă , atunci rezolvă .

Pasul 3. În caz contrar, fie o submulțime a mulțimii ce conține indicii a cel puțin unui negativ, și fie

Pasul 4. Dacă , atunci se înlocuiește cu , și se revine la Pasul 1.

Pasul 5. Dacă . atunci se înlocuiește cu . se înlocuiește cu și se revine la Pasul 1.

Nu este dificil de văzut că șirul valorilor obiectiv obținut prin această procedură este nedescrescător și că procedura se termină într-un număr finit de cicluri.

Evident, procedura de restrângere se găsește într-o corespondență duală cu relaxarea. Astfel nu este surprinzător faptul că ea este mai utilă în problemele cu multe variabile decât în cele cu multe restricții. Observăm că restricțiile vor conține deseori condițiile , acest lucru putându-se realiza întotdeauna prin transformări adecvate.

Atunci prin restrângerea se elimină pur și simplu variabilele din problemă. Bacă este liniară, și dacă se alege ca mulțime a indicilor variabilelor nebazice curente, atunci restrângerea devine identică cu metoda simplex primal (cel puțin dacă toate bazele primalei sunt nedegenerate).

În programele liniare cu multe coloane, restrângerea devine astfel o schemă de generare de coloane. Mulțimea din pasul 3 este atunci o submulțime de variabile nebazice având factori de cost relativ negativi, iar eficiența algoritmului depinde de cât de ușor poate fi găsită o astfel de submulțime. Desigur, pentru a reduce numărul de variabile nu este necesar ca să fie egală cu . Aplicarea restrângerii asupra unei submulțimi de p restricții liniar independente permite eliminarea a p variabile.

Deși procedeul de restrângere, la fel ca și relaxarea, și-a găsit cea mai largă aplicare la programele liniare, utilitatea lui nu se limitează numai la aceste probleme. Un exemplu excelent de aplicare a restrângerii la programe neliniare este dat de algoritmul de partiționare al lui Rosen din 7.3. Acolo restricțiile nu sunt liniare în toate variabilele, ceea ce conduce la o problemă neconvexă. De aceea, în afara pașilor prezentați aici, apar pași suplimentari, care includ o liniarizare a problemei .

2.2.2. Proceduri de partiționare în programarea neliniară

2.2.2.1. Introducere

La fel ea pentru programele liniare, termenul „procedură de partiționare" se va folosi aici pentru a desemna un algoritm care împarte variabilele problemei în două submulțimi, acționând la început asupra valorilor variabilelor dintr-una dintre submulțimi, apoi asupra celorlalte, variabile. Desigur, fiecare dintre submulțimi poate fi partiționată în continuare, cum ar fi cazul în care variabilele lor apar liniar și matricea lor restricție este bloc diagonală. Pentru problemele pe care le vom considera aici, partiționarea care se va folosi este în mod evident determinată de structura problemei.

Se vor discuta două clase de probleme de programare matematică înrudite. Prima clasă conține probleme de forma

să se minimizeze (2-74)

în raport cu (2-75)

Aici matricele , coeficienții de cost și vectorii din partea dreaptă sunt funcții continue de vectorul al variabilelor de cuplare. Pentru fixat problema este liniară în , astfel încât este natural ca variabilele și să fie ajustate separat. Probleme de acest tip reprezintă generalizarea pentru cazul neliniar a programelor liniare dual-unghiulare. Pentru a trata programe de forma (2-74) – (2-75), Rosen și-a extins procedeul de partiționare pentru programe liniare descris în 5.4. Această extindere va fi discutată în prima parte a capitolului.

A doua clasă de probleme considerata aici este de forma:

să se minimizeze (2-76)

în raport cu

(2-77)

(2-78)

unde și sunt funcții de cu valori solare, respectiv vectoriale și este o submulțime arbitrară a lui . Dată fiind libertatea în a alege pe , în această clasă de probleme pot fi incluse unele noi tipuri importante. Cel mai interesant caz apare atonei când se ia ca fiind mulțimea tuturor vectorilor – dimensionali cu componente întregi. Dacă, în plus, și sunt liniare, atunci (2-76) – (2-78) este binecunoscutul program liniar cu variabile mixte. Dacă este de formă bloc diagonală și , atunci (2-76) – (2-78) este un caz special al problemei (2-74) – (2-75), în care și nu sunt funcții de .

La fel ca în metoda lui Rosen, la început se modifică , apoi se rezolvă un program liniar pentru noi valori ale lui etc. Dar similitudinea între cele două proceduri se termină aici. Algoritmul lui Benders ajustează variabilele rezolvând un program pe , la care se adaugă la fiecare ciclu una sau mai multe restricții. Aceste restricții se construiesc folosind, pentru fixat, soluția duală a lui (2-76) – (2-78). Atunci, dacă (2-76) – (2-78) este un program liniar cu variabile mixte, programul care determină "variația lui este un program liniar, în întregime cu "variabile întregi. Astfel, în acest caz, (2-76) – (2-78) se rezolva modificând separat variabilele discrete și continue. Acest fapt prezintă interes deosebit atunci când A are structură specială, de exemplu este bloc diagonală sau de tip transport.

Deoarece acest program nu este în general convex, i se pot numai adăuga restricții, fără a fi posibilă eliminarea unora dintre ele. Metoda lui Rosen folosește restrângerea atunci când cere ca o submulțime de inegalități liniare cu matricea coeficienților nesingulară să fie satisfăcută cu egalitate. Aceasta este apoi folosită pentru a elimină variabilele în funcție de , ceea ce permite construirea unei probleme de coordonare în .

2.2.2.2. Algoritmul de partiționare al lui ROSEN pentru programe neliniare

Algoritmul de partiționare al lui Rosen pentru programe liniare, discutat în 5.4, a fost extins la o clasă de programe neliniare. Procedura este aplicabilă la problemele de forma

să se minimizeze (2-79)

în raport cu

(2-80)

Aici variabilele problemei se partiționează în două submulțimi: variabile de cuplare și variabile „liniare" . Vectorii coeficienților de cost , coeficienții matricei , și vectorii din membrul al doilea pot fi toți funcții liniare sau neliniare de variabilele .

O mare varietate de probleme poate avea această structură. Caracteristica deosebită a unei astfel de structuri constă în faptul că pentru fixat, problema se descompune în subprobleme liniare independente, în variabilele a. Astfel, o procedură iterativă care la început să modifice pe și apoi pe apare ca fiind naturală.

Matricele de bază și prețurile umbră ale soluțiilor subproblemelor se folosesc pentru formarea unui program coordonator în variabilele . Aici, spre deosebire de algoritmii binari, programul coordonator are o funcție obiectiv neliniară.

Convergența către un minim global în (2-79) – (2-80) se poate demonstra doar în cazul în care coeficienții și matricele sunt independente de , astfel încât problema devine

să se minimizeze (2-81)

în raport cu

(2-82)

Aceasta este forma pentru care vor fi formulate și demonstrate proprietățile algoritmului,urmând ca o forma mai generală a lui (2-79) – (2-80) să fie considerată ulterior. Se fac următoarele ipoteze suplimentare;

PREZENTAREA ALGORITMULUI

Problema (2-81) – (2-82) poate fi rezolvată minimizând la început în raport cu pentru toți (fixați), iar apoi minimizând rezultatul în raport cu . Putem astfel scrie

să se minimizeze (2-83)

unde este mulțimea tuturor valorilor lui pentru care fiecare minimizare din interiorul

parantezei este admisibilă. Dacă

(2-84)

Atunci (2-85)

Vectorii se vor numi admisibili. Se definește

(2-86)

Astfel problema (2-83) devine:

să se minimizeze (2-87)

în raport cu (2-88)

Următoarea teoremă afirmă că, în condițiile ipotezei 1, (2-87) – (2-88) este un program convex.

Dacă fiecare este convex, atunci mulțimea este convexă și fiecare funcție este convexă pe .Astfel problema originală a fost redusă la rezolvarea unui program convex în . Totuși, acest program apare ca fiind dificil datorită faptului că simpla evaluare a uneia dintre funcțiile necesită rezolvarea subproblemei liniare

să se minimizeze subproblema (2-89)

în raport cu (2-90)

Totodată, restricțiile asupra lui nu apar sub forma explicită a unor funcții. Pentru a atenua aceste dificultăți, vom încerca să obținem cât mai multe informații posibile din rezolvarea programului liniar (2-89) – (2-90). Aceste informații vor determina funcția și limitările lui în vecinătatea vectorului curent .

Fie vectorul admisibil și să presupunem că fiecare dintre programele liniare (2-89) – (2-90) a fost rezolvat pentru . Această rezolvare partiționează pe (care are mai multe linii decât coloane) într-o matrice de bază pătrată și nesingulară și într-o parte nebazică . Matricea determină vectorul soluție prin relația

(2-91)

Dacă admitem (temporar) că (2-91) determină relația între x și y, subproblema (2-89) – (2-90) poate fi scrisă astfel:

să se minimizeze (2-92)

în raport cu (2-93)

(2-93)

Relația (2-93) poate fi inversată pentru a exprima pe în raport cu :

(2-94)

Ecuația (2-94) este folosită pentru a elimina pe y din (2-92) și (2-93) obținând

(2-95)

(2-96)

Unde (2-97)

este vectorul variabilelor duale sau multiplicatorii simplex asociați bazei optimale .

Presupunând că are loc relația (2-94), va fi dat de (2-95), iar restricțiile (2-96) determină mulțimea din (2-84). Astfel în spațiul variabilelor poate fi enunțată o problemă de coordonare cu funcția obiectiv

(2-98)

și restricțiile (2-96). Această abordare ar putea fi considerată ca o aplicație a procedurii de restrângere schițată în 5.2, în care se cerea ca anumite restricții ale problemei originale [cele din (2-93)] să fie realizate cu semnul egal. Funcția obiectiv (2-98) este convexă, deoarece funcțiile sunt convexe și . Cu toate acestea, restricțiile (2-96) sunt neliniare și pot să nu determine o mulțime convexă, deoarece elementele lui pot avea orice semn.

Rezolvarea programelor cu restricții neliniare este în mod considerabil mai dificilă decât în cazul în care ele ar fi liniare, iar soluția globală nu poate fi garantată atât timp cât nu este asigurată convexitatea. Totuși, orice program coordonator trebuie doar să poată furniza un vector admisibil îmbunătățit , bineînțeles în cazul în care un astfel de vector există.

Având în vedere acest lucru, vom liniariza restricțiile (2-96) în vecinătatea punctului curent . Proprietățile programului coordonator care va rezulta sunt tratate în teorema 3. Fie

(2-99)

și definim matricea jacobiană corespunzătoare lui . Aceasta este matricea ale cărei linii sunt gradienții componentelor lui evaluați în punctul .

Liniarizând(2-96) în vecinătatea lui , se obține

(2-100)

Problema de coordonare are funcția obiectiv convexă (2-98) și restricțiile liniare (2-100)

să se minimizeze

problemă de coordonare (2-101)

în raport cu

(2-102)

Fie o soluție a acestei probleme. Deoarece este admisibil pentru problema de mai sus, vom avea

(2-103)

Desigur, (2-103) poate fi realizată cu egalitate.

În acest caz soluția curentă

(2-104)

poate fi testată pentru optimalitate.

Rezumat. Rezultă următorul rezumat al algoritmului:

Pasul 1. Se alege un vector .

Pasul 2. Se rezolvă subproblemele liniare (2-89) – (2-90) cu ,obținând soluții optimale matrice de bază și variabile duale

Pasul 3. Folosind cantitățile obținute la pasul 2, se formează problema de coordonare (2-101) – (2-102) care se rezolvă, obținând o soluție optimală .

Pasul 4. Dacă și toți din (2-104) sunt nenegativi, atunci soluția este optimală.

Pasul 5. Dacă și pentru , atunci, pentru subproblema liniară există o variantă de bază optimală în care cel puțin o linie a matricei corespunzătoare unei componente negative a lui este absentă. Se execută cel puțin o schimbare de bază și se revine la Pasul 3 cu noi matrice de bază, noi variabile duale etc.

Pasul 6. Dacă și toți atunci rezultă că – Se revine la pasul 2 cu înlocuit de .

Pasul 7. Dacă și un anumit , se calculează

a) Dacă , atunci se obține că .

Se revine la Pasul 2 cu înlocuit prin

b) Dacă , atunci, pentru anumiți și avem și .

Se poate arăta că în acest caz linia a matricei se poate schimba reciproc cu o anumită linie a lui furnizând o nouă bază optimală pentru subproblema liniară . Folosind această nouă bază, se revine la pasul 3.

Deoarece subproblemele liniare sunt scrise cu mai multe inegalități decât necunoscute, cel mai eficient pentru rezolvarea lor va fi un algoritm dual. Principala sa dificultate constă în mod obișnuit în numărul mare de restricții.

Metoda lui Rosen de proiecție a gradientului este astfel adecvată rezolvării problemei de coordonare, deoarece aceasta este o metodă duală a cărei eficiență depinde în principal de numărul de variabile al problemei și numai în mod secundar de numărul de restricții. Dacă numărul de iterații tinde către infinit, proiecția gradientului aplicată problemelor cu funcție obiectiv convexă și restricții liniare converge către soluția globală.

În ceea ce privește calculul, un aspect important al algoritmului este acela că fiecare iterație produce o soluție admisibilă pentru problema originală, cu valoare obiectiv necrescătoare. Deci calculele pot fi terminate la orice punct cu o soluție a cărei valoare obiectiv nu este mai mare decât toate cele anterioare.

2.3. Programe neliniare cu variabile de cuplare (caz particular)

2.3.1. Prezentare

O clasă importantă de programe neliniare are forma:

să se minimizeze

în raport cu (2-53)

Aici coeficienții de cost , matricele restricțiilor și vectorii din membrul al doilea sunt funcții de vectorul variabilelor de cuplare, . Dacă se dau valori acestor variabile de cuplare, atunci problema este liniară în vectorii și are o matrice a restricțiilor de formă bloc diagonală, putându-se astfel descompune în p subprobleme independente.

Problemele de acest tip sunt generalizări neliniare ale programelor liniare cu structură dual unghiulară, reducându-se la acest caz atunci când și sunt independente de y, iar și sunt liniare. Astfel de probleme se pretează la o rezolvare cu ajutorul metodelor de separare, în cadrul cărora se face la început ca y să varieze, apoi se rezolvă subproblemele liniare pentru a obține noi .

Astfel de probleme apar în mod frecvent în planificarea producției și distribuției produselor din industria petrolieră, unde se urmărește optimizarea globală ă unui complex de rafinării.

Fiecare rafinărie conține un număr de unități de proces (distilare, cracare, reformare etc.) care transformă intrările de surse brute în produse finite (benzină, motorină, ulei).

Acestea se trimit apoi la stațiile terminale de expediție, pentru a fi puse la dispoziția beneficiarilor. Trebuie să fie luate următoarele decizii tipice:

Cum trebuie alocate diferitele materii prime la fiecare rafinărie precum și la diferitele unități de proces din interiorul fiecărei rafinării ?

Cum se vor aloca fluxurile din cadrul fiecărei rafinării către unitățile de proces? Această decizie este complicată de faptul că anumite fluxuri intermediare dintr-o rafinărie pot fi trimise către alte rafinării în vederea unor utilizări viitoare.

Ce condiții de operare trebuie menținute în diferitele unități de proces? Acestea includ temperatura, presiunea, ratele de reciclare și activitățile de catalizare. Modul în care proprietățile fluxului depind de aceste variabile este aproape întotdeauna neliniar.

Ce materiale trebuie amestecate în fiecare produs pentru ca să se asigure proprietățile sortimentale corespunzătoare. Un exemplu de acest fel este necesitatea de a asigura o cifră octanică minimă pentru benzinele de motor. Și aici, componențele care trebuie amestecate pot fi transportate între rafinării.

Care dintre cererile de produse finite trebuie satisfăcute de fiecare rafinăriei

Aceste decizii sunt în general corelate. Deciziile de transport din 5 depind nu numai de costurile de transport dar și de disponibilitățile din diferitele produse la fiecare rafinărie 4 și de eficientele relative cu care pot fi produse diferitele produse 1—3. Deoarece rafinăriile diferă deseori semnificativ în ceea ce privește vârsta, capacitatea, echipamentul etc., problema poate fi destul de complexă. Perioada de timp luată în considerare se presupune a.fi destul de lungă, adică un sfert sau o jumătate de an, astfel încât complicațiile implicate de operațiile zilnice pot fi neglijate.

Figura 2 – 5 Complex de rafinării

2.3.2. Exemplu de abordare a unei probleme neliniare cu variabile de cuplare

Pentru a vedea cum o astfel de problemă conduce la un program neliniar, să considerăm următorul exemplu simplificat.

Două rafinării și, fabrică produsele și, pentru a fi livrate la stațiile terminale și . Se presupun cunoscute cererile cumulate pentru aceste produse la fiecare stație terminală de-a lungul intervalului de planificare; cererea de , la se notează . La rafinăria dispunem de trei materiale , și care amestecate dau și . În plus, poate fi transportat de la la unde se amestecă cu și .

Sistemul este reprezentat schematic în Figura 2 – 6. Să considerăm la început operația de amestec de la rafinăria . Produsele și trebuie să satisfacă cifrele octanice minime specificate și . Cifra octanică depinde atât de volumul fiecărui material prezent în amestec cât și (în mod neliniar) de concentrația tetraetilului de plumb (TEP) adăugat, . Fie cifra octanică a materialului atunci când concentrația TEP este y și fie volumul materialului folosit în amestec pentru a produce produsul .

În mod obișnuit se presupune contribuția lui la cifra octanică a produsului j este produsul lui cu volumul fracționar al lui , și contribuțiile diferitelor materiale sunt aditive.

Astfel, dacă este cifra octanică minimă a produsului , vom obține restricția:

(2-54)

care poate fi rescrisă sub forma:

. (2-55)

Figura 2 – 6 Sistem de proces și distribuție

Variabilele de amestec sunt supuse unei restricții prin care se exprimă faptul că nu se poate folosi dintr-un material mai mult decât este disponibil. Cantitățile disponibile sunt determinate de către unitățile de proces care întrețin operația de amestec, adică de către variabilele de control Dacă reprezintă disponibilitatea din materialul , atunci restricția de a nu se folosi din fiecare material mai mult decât este disponibil se scrie astfel:

(2-56)

În formă matriceală restricțiile vor fi scrise:

(tabelul 2 – 1)

unde . (2-57)

O mulțime similară de restricții corespunde amestecului din rafinăria 2 :

(tabelul 2 – 2)

unde este concentrația de TEP folosită la . Se observă că operația de transport a lui de la la a introdus o variabilă comună în problemele de amestec ale acestor rafinării, care este tratată ca o variabilă de cuplare.

Dacă am avea de transportat trei sau patru materiale, atunci va fi de preferat să se combine aceste două blocuri într-unul singur, decât să se introducă variabile de cuplare suplimentare. Aceasta se datorează faptului că în algoritmul de partiționare al lui Rosen efortul de calcul crește odată cu numărul variabilelor de cuplare.

Restricțiile care descriu operația de transport exprimă faptul că cererile trebuie să fie satisfăcute și că nu se poate transporta mai mult decât este disponibil. Forma lor matriceală este:

(tabelul 2 – 3)

Pentru a obține o mulțime completă de restricții trebuie scrise ecuațiile de balanță materială care exprimă faptul că volumul de trimis de la la terminale este egal cu volumul total de materiale, din care prin amestec a rezultat . Astfel, de exemplu,

(2-58)

Aceste restricții servesc la o nouă cuplare a sistemului — variabilele sunt variabile de cuplare suplimentare. În cazul în care ele sunt incluse în restricțiile de amestec ale rafinăriilor 1 și 2, sistemul general apare ca în (tabelul 2 – 1) Acesta are forma din (2-53) cu cinci variabile de cuplare și trei blocuri diagonale, fiecare reprezentând o fază distinctă a funcționării sistemului.

Funcția obiectiv a sistemului este venitul minus costul. Venitul este fix în cazul în care toate cererile trebuie satisfăcute exact, astfel încât trebuie minimizat costul. Costurile se referă la transport,

Transport de gazolină de la rafinării la terminale amestec și prelucrare. Dacă este costul unitar de transport de la , atunci costul total de transport este

, (2-59)

unde reprezintă o variabilă corespunzătoare unei celule din tabloul de transport (tabelul 3). Dacă este costul unității de volum de TEP, atunci drept cost al amestecului poate fi adoptat costul cantității totale de TEP folosită:

(2-60)

astfel încât coeficientul de cost al variabilelor de amestec depinde de . Putem reuni toate costurile de prelucrare, definind funcțiile , care reprezintă costul producerii unui volum din materialul , atunci când variabila de control fixată este Astfel, costul total de prelucrare devine

(2-61)

Costul total se obține adunând cu

Un model realist ar cere o reprezentare mult mai detaliata tuturor operațiilor în, particular a prelucrărilor anterioare amestecului. Probleme de această formă cu până la 60 variabile de cuplare, șase blocuri și 500—600 restricții inegalitate au fost rezolvate în industria petrolieră.

3. Metode de căutare evolutive

Dintre aceste metode voi exemplifica în continuare o subcategorie și anume metode de căutare evolutive. Algoritmii al căror mod de funcționare îi voi descrie sunt următorii:

Metoda NELDER – MEAD

Metoda BOX

Algoritmul PSO (Particle Swarm Optimisation)

3.1. Metode de căutare directă

Metodele de căutare directă pot oferi numai soluții aproximative ale problemelor de optimizare, dar în schimb asigură convergența pentru orice punct ales inițial, în timp ce metodele care folosesc gradientul, matricea Hessian sau o aproximație a acesteia – asemenea metode sunt denumite uneori metode indirecte – pot oferi soluții exacte, dar necesită un punct inițial bine ales.

În cazul metodelor indirecte, folosirea indicațiilor date de gradient (sau de gradient și de matricea Hessian), de exemplu pentru determinarea unui maxim, poate fi comparată cu o ascensiune spre un vârf de munte care este vizibil în timpul urcușului, iar folosirea metodelor directe poate fi considerată analoagă cu o ascensiune în condiții de lipsă de vizibilitate a vârfului (de exemplu, pe ceață), când singurele informații de care dispune cel care caută vârful se referă la sesizarea faptului că urcă sau coboară.

Figura 3 – 1

Ilustrarea metodelor directe poate fi ușor făcută prin intermediul unei funcții de o singură variabilă , variabila fiind deci scalară; o asemenea căutare a extremului este denumită „de-a lungul unei drepte", deoarece variabila ia valori în , deci pe o dreaptă.

Presupunând că se cunoaște intervalul de variație (al mărimii ) în care se găsește extremul – de exemplu se știe că între 0 și 1 există un maxim, funcția fiind concavă – căutarea directă poate fi efectuată prin mai multe metode. În Figura 3 – 1 este ilustrată aplicarea metodei dihotomiei (înjumătățirii intervalului).

Metoda prevede alegerea a două valori și în jurul centrului intervalului – respectiv în jurul valorii x = – cele două valori fiind apropiate între ele, deci diferind cu o valoare mică ; se compară valorile funcției și și dacă rezultă, de exemplu

(3-1)

ca în Figura 3 – 1, atunci poate fi eliminat intervalul (hașurat în Figura 3 – 1), întrucât funcția este concavă și maximul se va găsi între și .

Se aleg apoi alte două valori și în jurul centrului intervalului rămas (dintre și .0 ), separate tot prin , și se compară valorile și ale funcției criteriu; dacă rezultă, de exemplu,

(3-2)

ca în Figura 3 – 1, atunci poate fi eliminat și intervalul (hașurat) întrucât maximul se va găsi în domeniul respectiv în intervalul rămas nehașurat

(3-3)

Se continuă succesiv cu centrări de perechi de valori , comparații ale valorilor funcției și înjumătățiri de interval, până la obținerea unui interval foarte mic în care se găsește punctul de optim , deci până la obținerea soluției aproximative cu o toleranță admisă.

Neglijând valoarea a distanței dintre perechile de puncte și având în vedere că prin două testări ale valorilor funcției criteriu se asigură o înjumătățire a fiecărui interval, rezultă că după testări intervalul care rămâne de explorat are expresia aproximativă

. (3-4)

Alte metode de căutare de-a lungul unei drepte (de exemplu, metoda Fibonacci și metoda secțiunii de aur) au o eficiență sensibil superioară, reducerea intervalului inițial fiind mai rapidă.

Metodele de căutare de-a lungul unei drepte au o importanță deosebită și în căutarea extremului unei funcții criteriu de variabilă vectorială, întrucât după alegerea direcției din cadrul procesului iterativ, determinarea valorii se efectuează prin găsirea extremului unei funcții de variabila scalară a de-a lungul direcției .

Încercări de a organiza căutarea directă a punctului de optim, fără determinarea gradientului, au fost făcute și în cazul funcțiilor criteriu de variabilă vectorială, fiind elaborate diferite strategii de explorare. În cadrul celei mai simple strategii fiecare variabilă (componentă a vectorului ) este succesiv modificată până se obține un extrem – de exemplu un maxim – al funcției criteriu pe direcția variabilei respective, trecându-se apoi la modificarea variabilei următoare; procesul este oprit când nu se mai obțin creșteri ale funcției criteriu pe direcțiile nici unei variabile.

O strategie elaborată de Hooke, Jeeve și Wood prevede executarea dintr-un punct a unui „pas de explorare", reprezentând modificarea variabilelor cu o mică variație pozitivă prestabilită; dacă se constată o apropiere de optimul funcției criteriu, de exemplu de un maxim, atunci punctul obținut devine o „nouă bază" – spre deosebire de punctul de plecare reprezentând „vechea bază" – din care se efectuează pași de lucru pe direcțiile tuturor variabilelor (constituind un pas n-dimensional), pe fiecare direcție pașii fiind proporționali cu diferențele valorilor variabilei respective în noua și vechea bază.

Dacă pasul n-dimensional de lucru asigură apropierea de maxim, urmează un nou pas de explorare ș.a.m.d. Dacă pasul de lucru nu asigură creșterea funcției criteriu are loc anularea lui și executarea unui nou pas de explorare.

Dacă în direcția unei variabile pasul inițial de explorare nu determină o apropiere de maxim, atunci se efectuează un pas negativ de explorare pe aceeași direcție; în cazul când se obține o apropiere de maxim, acest pas rămâne valabil, dar dacă nici pasul negativ nu marchează creșterea funcției criteriu, atunci pe direcția respectivă nu se mai execută nici un pas de explorare.

Strategia elaborată de Nelder și Mead a fost numită „metoda simplex", într-un spațiu – dimensional un număr de puncte formând un „simplex". Pentru ilustrare, se consideră cazul minimizării unei funcții criteriu în care variabila vectorială are două componente, și , deci dependența poate fi reprezentată într-un spațiu . Întrucât asemenea reprezentări sunt dificile, se preferă intersecția corpului care reprezintă dependența cu o serie de plane paralele cu planul , prin intersecție rezultând în fiecare plan câte o curbă de valori constante ale funcției criteriu, numită și „curbă de nivel".

Proiectând diferitele curbe de nivel în planul variabilelor se obține o imagine de tipul celei din Figura 3 – 2, în care punctul de coordonate corespund minimului, deci

(3-5)

iar valorile constante – ale funcției criteriu pe diferitele curbe de nivel – sunt din ce în ce mai mici pe măsura apropierii de punctul , deci

(3-6)

Întrucât în un „simplex" se realizează cu puncte, în planul din Figura 3 – 2 un simplex va fi realizat cu 3 puncte, având forma unui triunghi. Se pornește de la un simplex inițial, de exemplu din Figura 3 – 2 și apoi se caută înlocuirea punctului cu cea mai mare valoare a funcției criteriu – deci cel mai depărtat de minim – de exemplu punctul , cu un alt punct de pe direcția , unde se găsește la jumătatea distanței dintre ; în. cazul unui număr mai mare de dimensiuni, punctul este centroidul tuturor vârfurilor rămase ale simplexului, după eliminarea vârfului .

Stabilind, printr-un algoritm adecvat de verificare a valorilor funcției criteriu, punctul de pe direcția care urmează să înlocuiască punctul – de exemplu punctul – se formează un nou simplex cu punctele continuându-se operațiile de construire a unei noi direcții pentru înlocuirea punctului cel mai depărtat de minim (din cadrul noului simplex) printr-un punct de pe noua direcție.

Figura 3 – 2

Când ajunge în apropierea minimului , simplexul are aspectul definit de punctele , cu direcția pentru înlocuirea punctului printr-un nou punct; în Figura 3 – 2, punctul coincide cu , ceea ce evident nu este obligatoriu.

Unele metode de căutare directă nu realizează o explorare prin stabilirea unui algoritm de elaborare a traseului, ca în cazurile menționate anterior, ci selectează în mod aleator punctele în care se determină valoarea funcției criteriu. Dacă punctele respective acopere complet gamele de variație estimate pentru toate variabilele de care depinde funcția criteriu, asemenea metode de căutare directă cu selectare aleatoare a punctelor de explorare sunt încadrate în clasa metodelor denumite „Monte Carlo".

În cadrul acestor metode, punctul în care a rezultat valoarea extremă a funcției criteriu este considerat ca optim. De fapt, după o primă explorare se poate stabili dacă optimul este localizat într-o anumită zonă restrânsă din cadrul celei considerate inițial și se poate calcula probabilitatea apropierii de optim, urmând apoi succesiv noi explorări aleatoare alternând cu noi reduceri ale zonei.

Elemente de calcul statistic intervin și în alte aspecte ale rezolvării problemelor de optimizare, nu numai în selectarea aleatoare a punctelor explorate.

Astfel, în unele cazuri însăși alegerea direcțiilor de căutare are un caracter aleator, rezultând metode aleatoare de căutare – prin procedee iterative – a optimului.

În alte cazuri, însăși variabilele de care depinde funcția criteriu (sau unii parametri care intervin în această funcție) sunt variabile aleatoare. în asemenea cazuri, optimizarea urmărește găsirea unui extrem al valorii medii a funcției criteriu, iar metodele de calcul folosite în acest scop sunt uneori denumite metode de programare stocastică.

3.2. Metoda NELDER—MEAD

3.2.1. Considerente teoretice

Este utilizată pentru optimizarea sistemelor multivariabile neliniare fără restricții și minimizează .

Esența metodei Nelder-Mead constă în deplasarea spre optim a figurii “N+1” dimensionale numită simplex, unde este dimensiunea spațiului variabilelor.

Deplasarea este asigurată prin operații de reflexie, expansiune sau contracție a simplexului pe baza testării valorilor funcției în vârfurile acestuia și a determinării celui mai favorabil și a celui mai defavorabil vârf.

Dacă prin testare se determină cele mai defavorabile vârfuri acestea sunt înlocuite prin imaginile lor reflectate care pot da succese sau insuccese pentru mișcarea spre optimul lui .

În cazul unor succese se procedează la o expansiune a simplexului îndreptată spre optimul funcției, iar la insuccese se realizează o contracție parțiala.

Când se ajunge în regiunea optimului simplexul se contractă în totalitate. Simplexul „N+1" dimensional este determinat de punctele care sunt vârfurile unui poliedru regulat din . Un simplex se generează alegând un punct de bază și o latură „a" cu care se construiesc celelalte vârfuri. Procedeul de generare a simplexului se bazează pe următoarele observații:

Distanța dintre oricare două puncte și , vârfuri alăturate ale simplexului este dată de relația:

(3-7)

Coordonatele unui vârf oarecare se pot scrie cu ajutorul punctului de bază inițial astfel:

(3-8)

Precizarea coordonatelor ale punctului se face cu ajutorul mărimilor intermediare și de valori:

(3-9)

(3-10)

Determinarea valorilor lui și se face presupunând că vârfurile fac parte din simplexul dimensional adică:

(3-11)

(3-12)

Cele vârfuri ale simplexului pot fi precizate cu următoarea reprezentare :

Centroidul simplexului față de care se realizează operațiile de reflexie, expansiune sau contracție se determină cu formula:

(3-13)

Printr-o strategie corespunzătoare simplexul definit de cele N + 1 vârfuri se „mișcă" spre punctul de optim al funcției .

Exemplu:

În un simplex este un triunghi iar în este un tetraedru

Figura 3 – 3

Figura 3 – 4

Reflectarea unui punct fata de centroid

Figura 3 – 5

Contracție (shrink) care apare foarte rar

Figura 3 – 6

Expansiunea unui punct față de centroid:

Figura 3 – 7

Contracția unui punct fata de centroid (în interior sau în exterior)

Figura 3 – 8

Figura 3 – 9

În realizarea simplexului este mai complicata dar principiul este același

Figura 3 – 10

Figura 3 – 11

Figura 3 – 12

Algoritmul metodei Nelder-Mead este următorul:

Pasul 1. Se alege un punct inițial care constituie vârful de bază al simplexului și latura „a".

Pasul 2. Se generează un simplex inițial prin vârfurile:

(3-14)

cu determinat prin valorile distribuite în tabloul prezentat mai sus.

Pasul 3. Se evaluează funcția criteriu în vârfurile simplexului. Cel mai defavorabil punct (în care funcția are valoarea cea mai mare) se înlocuiește cu un altul obținut prin una din operațiile de reflexie, contracție sau expansiune.

Pasul 4. Se realizează reflexia prin relația:

(3-15)

unde este o constantă pozitivă.

sunt coordonatele centroidului obținute astfel:

(3-16)

Pasul 5. Dacă punctul reflectat dă cea mai defavorabilă valoare a funcției criteriu se efectuează contracția prin relația:

(3-17)

unde 0<

Daca punctul reflectat este mai bun ca cel mai defavorabil punct dar nu este cel mai bun, se efectuează următoarea contracție:

(3-18)

Se evaluează funcția în acest punct și dacă se obține o îmbunătățire față de punctele curente se reia procedura cu următorul cel mai rău punct.

Dacă nu este obținută o îmbunătățire, vârfurile centroidului sunt contractate spre cel mai bun punct cu jumătatea distanței de la fiecare vârf la cel mai bun punct după relația:

(3-19)

Pasul 6. Dacă punctul reflectat calculat la pasul 4 este cel mai bun punct, se efectuează o expansiune după relația:

(3-20)

unde este o constantă pozitivă.

Dacă punctul de expansiune este o îmbunătățire față de punctul reflectat, se reține punctul de expansiune și procesul este reluat.

Dacă punctul de expansiune nu este o îmbunătățire față de punctul reflectat se reține punctul reflectat și se începe un nou stadiu.

Pasul 7. Căutarea încetează când este satisfăcut următorul criteriu de convergență:

(3-21)

unde

= vârfurile simplexului în afara celui mai defavorabil

= centroidul simplexului

Un exemplu de evoluție a simplexului pentru o funcție data în Figura 3 – 13

Figura 3 – 13

Organigrama algoritmului este prezentată în Figura 3 – 14.

Figura 3 – 14 Organigrama algoritmului Nelder –

3.2.2. Aplicație la metoda Neder – Mead

Figura 3 – 15

Figura 3 – 16

Figura 3 – 17

Figura 3 – 18

Figura 3 – 19

Figura 3 – 20

Figura 3 – 21

3.3. Metoda BOX

3.3.1. Considerente teoretice

Este utilizată pentru optimizarea sistemelor multivariabile, neliniare cu restricții pentru minimizarea funcției cu restricțiile:

(3-22)

Variabilele implicite sunt funcții de variabilele explicite .

Metoda este o extindere a metodei simplex (Nelder-Mead) pentru cazul optimului cu restricții al funcțiilor multivariabile.

Această metodă utilizează un simplex flexibil numit „complex" cu un număr de vârfuri. Complexul își modifică forma astfel încât să se încadreze în domeniul admisibil al restricțiilor și să înainteze spre optim.

În fapt se construiește un spațiu rectangular prin relațiile care dau restricțiile explicite de forma:

(3-23)

în care se desfășoară ajustări ale complexului pentru aducerea vârfurilor în domeniul admisibil și apropierii de optim.

Complexul se construiește dintr-un punct inițial de pornire admisibil, cele vârfuri suplimentare generându-se după relația:

(3-24)

unde este un șir de numere uniform distribuite sau aleatoare din intervalul (0,1).

Vârfurile complexului trebuie să satisfacă atât restricțiile explicite cât și pe cele implicite. Dacă restricțiile explicite sunt violate punctul este deplasat până intră în domeniul rectangular admisibil.

Dacă restricțiile implicite sunt violate de către un anumit vârf, punctul corespunzător acestui vârf se contractă cu jumătatea distanței până la centroidul complexului determinat de relațiile:

(3-25)

Rezultă că punctul contractat are coordonatele:

(3-26)

unde :

= coordonatele punctului contractat;

= coordonatele punctului – din afara domeniului restricțiilor.

Acest proces se repetă până ce toate vârfurile respectă restricțiile implicite.

Înaintarea spre punctul de optim se face evaluând funcția în vârfurile complexului și estimând vârfurile importante (cele mai defavorabile și cele mai favorabile) și funcție de aceste estimări se modifică complexul prin reflexii și contracții succesive.

Algoritmul Box este următorul:

Pasul 1. Se generează un complex inițial prin alegerea unui vârf și generarea celor vârfuri suplimentare cu relația:

(3-27)

și se calculează centroidul inițial cu relația:

(3-28)

Pasul 2. Se testează dacă vârfurile satisfac restricțiile explicite și dacă unele nu sunt satisfăcute de anumite vârfuri acestea se readuc în interiorul restricțiilor violate.

Pasul 3. Se testează dacă vârfurile satisfac restricțiile implicite și dacă unele nu sunt satisfăcute de anumite vârfuri acestea se reduc față de centroid după relația:

(3-29)

Pasul 4. Pentru complexul cu toate vârfurile în domeniul admisibil se evaluează în fiecare vârf.

Punctul în care funcția are valoarea cea mai defavorabilă este înlocuit cu un punct reflectat față de centroid care se determină prin relația:

(3-30)

S-a considerat = 1,3.

Pasul 5. Dacă punctul reflectat dă încă valoarea cea mai defavorabilă se realizează o contracție după o relație de tipul:

(3-31)

Pasul 6. Se cercetează dacă noul punct respectă restricțiile și se procedează ca mai înainte dacă restricțiile sunt violate.

Pasul 7. Se verifică criteriul de convergență exprimat prin:

pentru toate vârfurile și pentru un număr consecutiv de iterații:

Organigrama algoritmului este prezentată în Figura 3 – 22:

Figura 3 – 22

3.3.2. Aplicație la Metoda Box (Complex)

Figura 3 – 23

n=2 m=2 k=6

nr_max_iter=500

nr_restr__expl=2

alfa=1.300000

beta=0.000100

gama=5

delta=0.001000

––––––––––-

numere aleatoare:

r(1,1)=0.8462

r(1,2)=0.5252

r(2,1)=0.2026

r(2,2)=0.6721

r(3,1)=0.8381

r(3,2)=0.0196

r(4,1)=0.6813

r(4,2)=0.3795

r(5,1)=0.8318

r(5,2)=0.5028

r(6,1)=0.7095

r(6,2)=0.4289

coordonatele complexului inițial:

x[1][1]=16.9244

x[1][2]=10.5030

x[2][1]=4.0529

x[2][2]=13.4427

x[3][1]=16.7624

x[3][2]=0.3928

x[4][1]=13.6255

x[4][2]=7.5896

x[5][1]=16.6359

x[5][2]=10.0563

x[6][1]=14.1894

x[6][2]=8.5778

valorile funcției:

f[1]=0.2873

f[2]=-1.1079

f[3]=0.3510

f[4]=1.3938

f[5]=0.3161

f[6]=0.7636

––––- iterația 1–––-

coord. punctului corectat :

x[2][1]=17.8133

x[2][2]=3.7125

f[1] = 0.2873

f[2] = 1.5338

f[3] = 0.3510

f[4] = 1.3938

f[5] = 0.3161

f[6] = 0.7636

coordonatele centroidului :

x_centr[1] =15.6275

x_centr[2] =7.4239

–––– iterația 2 –––-

coord. punctului corectat :

x[1][1]=14.3504

x[1][2]=0.2974

f[1] = 0.4200

f[2] = 1.5338

f[3] = 0.3510

f[4] = 1.3938

f[5] = 0.3161

f[6] = 0.7636

coordonatele centroidului :

x_centr[1] =15.8053

x_centr[2] =6.0658

–––– iterația 3 –––

coord. punctului corectat :

x[5][1]=14.5112

x[5][2]=2.0575

f[1] = 0.4200

f[2] = 1.5338

f[3] = 0.3510

f[4] = 1.3938

f[5] = 1.6426

f[6] = 0.7636

coordonatele centroidului :

x_centr[1] =15.3482

x_centr[2] =4.1140

–––– iterația 4 –––-

coord. punctului corectat :

x[3][1]=13.6861

x[3][2]=7.0822

f[1] = 0.4200

f[2] = 1.5338

f[3] = 1.9629

f[4] = 1.3938

f[5] = 1.6426

f[6] = 0.7636

coordonatele centroidului :

x_centr[1] =14.8980

x_centr[2] =4.4470

–––– iterația 5 –––-

coord. punctului corectat :

x[1][1]=14.8999

x[1][2]=7.5936

f[1] = 1.9403

f[2] = 1.5338

f[3] = 1.9629

f[4] = 1.3938

f[5] = 1.6426

f[6] = 0.7636

coordonatele centroidului :

x_centr[1] =14.7651

x_centr[2] =5.8039

–––– iterația 65 –––

coord. punctului corectat :

x[4][1]=15.0053

x[4][2]=4.9842

f[1] = 3.9866

f[2] = 3.9866

f[3] = 3.9866

f[4] = 3.9866

f[5] = 3.9866

f[6] = 3.9866

coordonatele centroidului :

x_centr[1] =15.0144

x_centr[2] =4.9786

–––– iterația 66 –––

coord. punctului corectat:

x[6][1]=15.0209

x[6][2]=4.9935

f[1] = 3.9866

f[2] = 3.9866

f[3] = 3.9866

f[4] = 3.9866

f[5] = 3.9866

f[6] = 3.9866

coordonatele centroidului :

x_centr[1] =15.0138

x_centr[2] =4.9811

–––––––––––-

valoarea finala a funcției :

f[6]=3.9866

valorile finale ale lui x :

x[6][1]=15.0209

x[6][2]=4.9935

stop !

Figura 3 – 24

Ultimul pas:

–––– iterația 96 –––

coord. punctului corectat :

x[6][1]=15.0233

x[6][2]=4.9817

Valorile funcției în toate vârfurile:

f[1] = 3.9867

f[2] = 3.9867

f[3] = 3.9867

f[4] = 3.9866

f[5] = 3.9867

f[6] = 3.9867

coordonatele centroidului :

x_centr[1] =15.0194

x_centr[2] =4.9781

–––––––––––-

valoarea finala a funcției :

f[2]=3.9867

valorile finale ale lui x :

x[2][1]=15.0100

x[2][2]=4.9819

Figura 3 – 25

Ultimul pas:

–––– iterația 75 –––

coord. punctului corectat :

x[4][1]=15.0215

x[4][2]=4.9706

Valorile funcției în toate vârfurile:

f[1] = 3.9866

f[2] = 3.9866

f[3] = 3.9866

f[4] = 3.9866

f[5] = 3.9866

f[6] = 3.9866

coordonatele centroidului :

x_centr[1] =15.0176

x_centr[2] =4.9777

–––––––––––-

valoarea finala a funcției :

f[3]=3.9866

valorile finale ale lui x :

x[3][1]=15.0151

x[3][2]=4.9706

stop !

Figura 3 – 26

Ultimul pas:

–––– iterația 78 –––

coord. punctului corectat :

x[4][1]=15.0096

x[4][2]=4.9914

Valorile funcției în toate vârfurile:

f[1] = 3.9867

f[2] = 3.9867

f[3] = 3.9867

f[4] = 3.9867

f[5] = 3.9867

f[6] = 3.9867

coordonatele centroidului :

x_centr[1] =15.0085

x_centr[2] =4.9861

–––––––––––-

valoarea finala a funcției :

f[5]=3.9867

valorile finale ale lui x :

x[5][1]=15.0110

x[5][2]=4.9817

stop !

Figura 3 – 27

Ultimul pas:

–––– iterația 85 –––

coord. punctului corectat :

x[2][1]=15.0077

x[2][2]=4.9889

Valorile funcției în toate vârfurile:

f[1] = 3.9867

f[2] = 3.9867

f[3] = 3.9866

f[4] = 3.9866

f[5] = 3.9866

f[6] = 3.9866

coordonatele centroidului :

x_centr[1] =15.0111

x_centr[2] =4.9842

–––––––––––-

valoarea finala a funcției :

f[1]=3.9867

valorile finale ale lui x :

x[1][1]=15.0097

x[1][2]=4.9797

stop !

3.4. Algoritmului PSO (Particle Swarm Optimisation)

3.4.1. Considerente teoretice

Particle Swarm Optimization este un model introdus în 1995 de către James Kennedy și Russell Eberhart, inspirat din comportamentul social al stolurilor de păsări și al bancurilor de pești. Conceput inițial ca "O metodă de optimizare a funcțiilor continue neliniare" diferită față de algoritmii tradiționali de căutare, Particle Swarm Optimization lucrează cu o populație de potențiale soluții din spațiul de căutare. Prin cooperare și competiție între indivizi această tehnică poate da rezultate mult mai bune când este aplicată problemelor complexe de optimizare.

Principala dificultate întâlnită în tehnicile de optimizare multimodală și în special în Particle Swarm Optimization este convergența prematură care se concretizează în soluții sub-optimale. Această lucrare propune și analizează o serie de tehnici ce tratează problema convergenței premature și aduc de asemenea îmbunătățiri în ceea ce privește diminuarea costurilor algoritmului.

Tehnicile propuse prevăd modificări în modul de inițializare a populației, tratarea adaptivă a unor parametri, continuarea căutării după ce algoritmul converge prin reinițializare la scară redusă și formularea condițiilor de oprire. Spre deosebire de alte metode ce particularizează algoritmul pentru anumite probleme, cele propuse în acesta lucrare sunt aplicabile pe varianta de bază a algoritmului și aduc rezultate mai bune indiferent de problema pe care trebuie să o rezolve.

Considerând un grup cu elemente intr-un spațiu de căutare de dimensiune . La iterația a algoritmului, vectorul de poziție alcătuit din elemente pentru fiecare particula, se calculează de forma:

, (3-32)

unde reprezintă vectorul de viteza, obținut la fiecare iterație prin următoarea formula:

(3-33)

Coeficientul de inerție este coeficientul de inerție, iar este vectorul viteza de la iterația k. Acest vector viteza este alcătuit din suma vectoriala a 2 vectori (tendința cognitiva) și (tendința sociala). Semnificația notațiilor este următoarea: reprezintă vectorul poziției cea mai bune a particulei la iterația , în timp ce reprezintă vectorul poziției cea mai bune al tot grupului pana la iterația .

Cei 2 coeficienți r1 și r2 sunt 2 numere aleatoare, în intervalul [0, 1], care dau caracterul stocastic al algoritmului. Cei 2 coeficienți, și sunt constant reale aleatorii, numiți coeficienți de corecție, alese în general în domeniul [0, 2].

Alegând cele 2 constante egale (), contribuțiile celor 2 tendințe (cognitive și sociala) vor fi egale. Cei 2 vectori dau direcția și deplasarea elementului la pasul , folosind poziția cea mai buna dintre toate particulele, dar și poziția globala a grupului.

Parametrii algoritmului PSO

Viteza

Viteza intr-un algoritm este limitata pentru a stabiliza algoritmul.

La o iterație k, viteza nou calculata poate atinge valori foarte mari, putând perturba algoritmul, putând trece ușor peste punctul optim. Pentru a evita astfel de divergente și pentru a controla algoritmul, fiecare particula are viteza limitata la o viteza maxima: ().

Figura 3 – 28

Figura 3 – 29

Viteza maxima nu trebuie limitata la o valoare prea mica; în acest caz am explorarea domeniului s-ar realiza prea încet.

Există doua modalități de a impune restricții asupra vitezei: în primul caz (a) se pot impune restricții asupra componentelor vitezei pe cele doua axe din sistemul cartezian al domeniului de căutare, sau, intr-o alta varianta (b), se pot pune restricții asupra normei vectorului, precum este indicat în figurile de mai sus:

Inițializarea poziției particulelor

Conform algoritmului, poziția fiecărei particulele este inițializată aleator în domeniul de căutare, pentru a se acoperi cât mai bine spațial de căutare.

Realizarea acestei inițializări în tot domeniul de căutare acorda algoritmului o diversitate mărită, și o probabilitate de descoperire a minimului funcției crescuta.

Factorul de inerție

S-a propus modificare factorului de inerție în funcție de numărul iterației pentru fiecare particular, în funcție de valoarea funcției pe care o au în acel moment.

Dacă o particula are valoarea funcției cost mai buna de cât cea a unei particule , probabilitatea ca optimul global să se afle în vecinătatea particulei este mai mare. Deci aceasta nu are nevoie de o inerție mare. În acest caz, particula trebuie să exploateze mai cu atenție spațial în care se afla, să nu se îndepărteze, pe când particula trebuie să exploateze alte regiuni.

Astfel factorul de inerție este definit corespunzător:

Se calculează un factor index pentru fiecare particula i

(3-34)

unde și sunt pozițiile particulei cu cea mai proasta, respectiv cea mai buna valoare a funcției calculate la iterația pentru particula

Se calculează factorul de inerție pentru fiecare particula

(3-35)

unde și sunt factorii de inerție cel mai mic, respectiv cel mai mare la iterația pentru particula .

Acest calcul adaptiv al factorului de inerție acorda un grad mai mare de flexibilitate a particulelor mai puțin performante, mărindu-le viteza, le da un spațiu mai mare de căutare a punctului optim.

Scopul strategiilor evoluționiste este de a îmbunătăți procesul de căutare aleatorie intr-o problema de optimizare (inițializarea problemei). Prin asta se înțelege identificarea domeniului care conține soluția optima. Presupunem deci ca suntem intr-un domeniu care permite găsirea soluției problemei. Exista diverse strategii de identificare și explorare a domeniilor cu probabilitate ridicata de găsire a soluției. Aceste strategii se împart în algoritmi genetici, evoluționiști, etc.

Algoritmul PSO a fost prezentat de Kennedy și Ederhart ca un algoritm stocastic de optimizare, fără calculul gradientului. Principiul fundamental al acestui algoritm este folosirea în comun a informației intre particule: la început, un grup de particule este distribuit intr-un domeniu de dimensiune , pentru ca mai apoi, particulele își calculează singure poziția în domeniu de căutare folosind o formulă simplă.

Poziția fiecărei particule la fiecare iterație este evaluata utilizând funcția obiectiv sau funcția cost. Memoria cognitiva a fiecărei particule permite găsirea celei mai bune valori a funcției cost. Mișcarea particulelor se poate face intr-un domeniu continuu sau unul discret.

Formula originala a lui Kennedy și Eberhart este următoarea:

(3-36)

Este considerat cazul general a unei distribuții intr-un domeniu , iar reprezintă un număr aleatoriu. Constantele și reprezintă ponderile celor 2 factori: cognitive și social.

Numerele aleatoare pot fi numere reale care să acorde ponderi diferite celor 2 factori sau pot fi vectori care măsoară fiecare componenta, cognitiva și sociala, pentru toate particulele din grup. Se observa caracterul aleatoriu al comportamentului grupului.

Formularea matematica

În formularea algoritmului, vectorul stocastic este dat de formula:

(3-37)

unde și reprezintă doua numere reale, scalare, aleatorii în intervalul [0,1]. Acestea măsoară ponderea celor doi factori, cognitiv și social.

// Initialize the particle positions and their velocities

for I = 1 to number of particles n do

for J = 1 to number of dimensions m do

X[I][J] = lower limit + (upper limit – lower limit) * uniform random number

V[I][J] = 0

enddo

enddo

// Initialize the global and local fitness to the worst possible

fitness_gbest = inf;

for I = 1 to number of particles n do

fitness_lbest[I] = inf

enddo

// Loop until convergence, in this example a finite number of iterations chosen

for k = 1 to number of iterations t do

// evaluate the fitness of each particle

fitness_X = evaluate_fitness(X)

// Update the local bests and their fitness

for I = 1 to number of particles n do

if (fitness_X[I] < fitness_lbest[I])

fitness_lbest[I] = fitness_X[I]

for J = 1 to number of dimensions m do

X_lbest[I][J] = X[I][J]

enddo

endif

enddo

// Update the global best and its fitness

[min_fitness, min_fitness_index] = min(fitness_X)

if (min_fitness < fitness_gbest)

fitness_gbest = min_fitness

for J = 1 to number of dimensions m do

X_gbest[J] = X(min_fitness_index,J)

enddo

endif

// Update the particle velocity and position

for I = 1 to number of particles n do

for J = 1 to number of dimensions m do

R1 = uniform random number

R2 = uniform random number

V[I][J] = w*V[I][J]

+ C1*R1*(X_lbest[I][J] – X[I][J])

+ C2*R2*(X_gbest[J] – X[I][J])

X[I][J] = X[I][J] + V[I][J]

enddo

enddo

enddo

Într-o variantă simplificată, algoritmul PSO este descris în pseudo-cod astfel:

Pentru fiecare particula din grup

Inițializează particula

Pana când este atins nr. maxim de iterații sau criteriul de eroare, repeta

Pentru fiecare particula din grup

Calculează funcția cost

Daca funcția cost actuala este cea mai buna a particulei (pbest)

Păstrează poziția și reactualizează pbest

Alege particula cu cea mai bună funcție cost din grup (gbest)

Pentru fiecare particula din grup

Calculează viteza particulei conform ecuației:

Actualizează poziția actuala a particulei conform ecuației:

O alta maniera de a generaliza algoritmul este data de folosirea în ecuație a coeficienților. Numerele aleatoare din ecuația vitezei pot fi păstrate intr-o matrice diagonala, pentru a ușura calculul vitezei pentru tot vectorul viteza:

(3-38)

, 0<<1 (3-39)

Algoritmul poate fi folosit și intr-un spațiu cu mai multe dimensiuni. Un algoritm PSO pentru un grup format din particule, intr-un spațiu de căutare de dimensiuni, este descris de următorul pseudo-cod:

Pentru fiecare particula ( )

Pentru fiecare dimensiune ( j = 1:n )

Calculează viteza particulei conform ecuației:

Actualizează poziția curenta a particulei conform ecuației:

3.4.2. Aplicație la metoda PSO

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

inertia = 1; %coeficient de inerție

correction_factor = 2; %factor de corecție

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Alegerea parametrilor

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAZ 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

inertia = 1; %coeficient de inerție

correction_factor = 3; %factor de corecție

swarm_size = 49; %numărul de puncte

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Figura 3 – 30

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAZ 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

inertia = 1; %coeficient de inerție

correction_factor = 4; %factor de corecție

swarm_size = 49; %numărul de puncte

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Figura 3 – 31

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAZ 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

inertia = 1.5; %coeficient de inerție

correction_factor = 2; %factor de corecție

swarm_size = 49; %numărul de puncte

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Figura 3 – 32

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAZ 4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

inertia = 1.9; %coeficient de inerție

correction_factor = 2; %factor de corecție

swarm_size = 49; %numărul de puncte

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Figura 3 – 33

În continuare, vom studia aplicarea algoritmului pentru diverse funcții.

În implementarea algoritmului vitezele au fost inițializate cu 0, c1= c2= 2, s-a considerat un grup format din p=49 de particule, iar calculele s-au efectuat pentru diverși factori de inerție w; numărul de iterații maxime s-a considerat 50.

Funcție simplă, în 2 dimensiuni

Pentru studiul acestei funcții am considerat următorii parametri:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

inertia = 1; %coeficient de inerție

correction_factor = 2; %factor de corecție

swarm_size = 49; %numărul de puncte

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Se poate observa că particulele grupului se mișcă uniform către găsirea punctului optim.

Eficienta algoritmului este demonstrata în evoluția grupului în funcție de numărul iterației. După 5 iterații este obținut minimul funcției și după 17 iterații toate particulele ajung în minimul funcției

Alegerea parametrilor

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAZ 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

inertia = 1; %coeficient de inerție

correction_factor = 0.2; %factor de corecție

În acest caz, datorita factorului de corecție mic, se poate observa ca punctul optim nu este găsit particulele grupului. Grupul se mișca uniform, dar prea încet pentru a fi eficient.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAZ 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

inertia = 2; %coeficient de inerție

correction_factor = 1; %factor de corecție

Datorita inerției mari, punctul optim este atins rapid, dar particulele oscilează mult în jurul lui. Grupul are o mișcare haotica.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAZ 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

inertia = 1; %coeficient de inerție

correction_factor = 2; %factor de corecție

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Datorită factorului de corecție mai mare, punctul optim este atins rapid, atât de către o particular, cât și de întreg grupul.

Funcție simpla, în 3 dimensiuni

Evoluția grupului în etape

Pentru studiul acestei funcții am considerat următorii parametri:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

inertia = 1; %coeficient de inerție

correction_factor = 2; %factor de corecție

swarm_size = 343; %numărul de puncte

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Eficiența algoritmului este demonstrată în evoluția grupului în funcție de iterație. Se observă că după 4 iterații se obține minimul funcției.

3.4.3. Concluzii pentru algoritmul PSO

O inerție prea mare provoacă o mișcare haotica a particulelor ce compun grupul. Deși în cazul unei inerția mari se poate ca valoarea optima să fie găsită mai repede, particulele oscilează foarte mult fata de punctul de optim, neputând să-l atingă.

În urma experimentului un factor de inerție 1.0 duce la o mișcare a particulelor uniforma către valoarea optima.

Un factor de corecție mic (0.1) duce la o întârziere a găsirii soluției optime foarte mult. Pe când un factor de corecție mare (2.0) determina găsirea soluției optime rapid, dar nu toate particulele ating repede acest punct optim.

În urma experimentului un factor de corecție 2.0 duce la găsirea punctului optim cât mai rapid de către toate particulele grupului.

3.4.4. Aplicație: Folosire ToolBox Matlab PSO

Am folosit aplicația intitulată ToolBox Matlab PSO în scopul de a compara rezultatele obținute de acestă rutină cu cele din rutina implementată de mine

Static test functions, minima don't change w.r.t. time/iteration:

1) Ackley

2) Alpine

3) DeJong_f2

4) DeJong_f3

5) DeJong_f4

6) Foxhole

7) Griewank

8) NDparabola (for this demo N = 2)

9) Rastrigin

10) Rosenbrock

11) Functia mea de test

12) Schaffer f6 modified (5 f6 functions translated from each other)

13) Tripod

Dynamic test functions, minima/environment evolves over time/iteration:

14) f6_bubbles_dyn

15) f6_linear_dyn

16) f6_spiral_dyn

Choose test function ? 11

1) Intense graphics, shows error topology and surfing particles

2) Default PSO graphing, shows error trend and particle dynamics

3) no plot, only final output shown, fastest

Choose plotting function ? 1

0) Minimize

1) Maximize

Choose search goal ? 0

Max velocity divisor (2 is a good choice) ? 2

How many particles (24 – 30 is common)? 24

0) Common PSO – with inertia

1) Trelea model 1

2) Trelea model 2

3) Clerc Type 1" – with constriction

Choose PSO model ? 0

PSO: 1/400 iterations, GBest = -2.185993353950098e-008.

PSO: 2/400 iterations, GBest = -2.185993353950098e-008.

PSO: 3/400 iterations, GBest = -2.4055572762057409.

PSO: 4/400 iterations, GBest = -2.4055572762057409.

PSO: 5/400 iterations, GBest = -2.4055572762057409.

PSO: 6/400 iterations, GBest = -2.8262822784152095.

…………………………………………………………..

PSO: 21/400 iterations, GBest = -3.7670365469090417.

PSO: 22/400 iterations, GBest = -3.7672117454472702.

…………………………………………………………

PSO: 365/400 iterations, GBest = -3.9866961272053061.

PSO: 365/400 iterations, GBest = -3.9866961272053061.

–> Solution likely, GBest hasn't changed by at least 1e-099 for 100 epochs.

Best fit parameters:

cost = f6( [ input1, input2 ] )

–––––––––––

input1 = 15.0162

input2 = 4.9837

cost = -3.9867

mean cost = -3.9208

# of epochs = 365

Figura 3 – 34

Figura 3 – 35

Figura 3 – 36

Figura 3 – 37

Figura 3 – 38

function [out]=f6(in)

x=in(:,1);

y=in(:,2);

1) Intense graphics, shows error topology and surfing particles

2) Default PSO graphing, shows error trend and particle dynamics

3) no plot, only final output shown, fastest

Choose plotting function ? 2

Figura 3 – 39

Figura 3 – 40

Best fit parameters:

cost = f6( [ input1, input2 ] )

–––––––––––

input1 = 15.0162

input2 = 4.9837

cost = -3.9867

mean cost = -3.957

# of epochs = 326

1) Intense graphics, shows error topology and surfing particles

2) Default PSO graphing, shows error trend and particle dynamics

3) no plot, only final output shown, fastest

Choose plotting function ? 1

0) Minimize

1) Maximize

Choose search goal ?1

Max velocity divisor (2 is a good choice) ? 2

How many particles (24 – 30 is common)? 24

0) Common PSO – with inertia

1) Trelea model 1

2) Trelea model 2

3) Clerc Type 1" – with constriction

Choose PSO model ? 0

Figura 3 – 41

Figura 3 – 42

Figura 3 – 43

Figura 3 – 44

Best fit parameters:

cost = f6( [ input1, input2 ] )

–––––––––––

input1 = 14.9718

input2 = 20.0282

cost = 4.9642

mean cost = 4.8555

# of epochs = 400

>>

4. SISCON – produs software pentru optimizare

4.1. Descriere soft (Manual de utilizare)

I. Descrierea pachetului de programe SISCON

Pachetul de programe SISCON, este o aplicație destinată optimizării și conducerii eficiente a sistemelor, dezvoltată în C++ (Builder).

Pentru ca în program să poată fi introduse orice tip de funcții sau combinații ale acestora, s-a realizat un analizor sintactic ce preia funcțiile ca șiruri de caractere, înlocuiește variabilele x prin seturi de valori numerice și calculează valorile funcției criteriu.

II. Lansarea sistemului SISCON

Fereastra de START a programului este:

Pentru a trece mai departe se face click în dreptunghiul interior ce determină apariția unui meniu de navigare arborescent, având ca scop ghidarea utilizatorului pentru a ajunge la metodele dorite. Prima fereastră ce apare după cea de start, are forma:

III Subrutine de optimizare pentru evaluarea deciziilor

Metode de optimizare multivariabile

III.1 Metoda Nelder – Mead

Metoda Nelder – Mead rezolvă probleme de tipul: min F(x1,x2,..,xn), unde F este o funcție neliniara, fără restricții.

III.2 Metoda Box

Metoda Box rezolvă probleme de tipul: min F(x1,x2,..,xn), unde F este o funcție neliniară cu restricții impuse variabilelor.

4.2. Metode de căutare evolutive

4.2.1. Aplicație la Metoda Nelder-Mead

4.2.1.1. Exemplul 1:

Figura 4 – 1

Figura 4 – 2

4.2.1.2. Exemplul 2:

Figura 4 – 3

Figura 4 – 4

4.2.2. Aplicație la Metoda BOX

Exemplu

Figura 4 – 5

Figura 4 – 6

4.3. Transformări echivalente pentru probleme de mari dimensiuni

4.3.1. Tehnici de descompunere

4.3.2. Tehnici de penalizare

4.3.3. Transformarea duală

5. Studiu comparativ al algoritmilor evolutivi

5.1. Exemple

Exemplul I Fie funcția:

Figura 5 – 1

Exemplul al II-lea

Funcție simplă, în 2 dimensiuni

Figura 5 – 2

Exemplul al III-lea:

Funcție simplă, în 3 dimensiuni

5.2. Observații

Pentru Exemplul I pot fi trase următoarele observații

Nelder-Mead

% Cu următoarea inițializare, minimul global va fi in punctul (20,15)

% deoarece se pleacă dintr-un punct apropiat de un minim local

%P(:,:,1)=[15 20 15;

% 15 15 20];

În testarea algoritmului am folosit diverse condiții inițiale și s-a observat că pentru alegerea simplexului cu vârfurile V1(15,15), V2(20,15) și V3(15,20) minimul obținut este unul local și nu unul global.

Din acest lucru rezultă și se confirmă ipoteza că este foarte important simplexul inițial.

BOX

Și pentru metoda Box a fost găsit un punct (coordonatele centroidului inițial) pentru care vom ajunge în același minim local (20,15) determinat cu metoda Nelder-Mead iar inițializarea este x=1, y=1;

PSO (Particle Swarm Optimisation)

De fiecare data s-a obținut soluția optimă.

Acest lucru se datorează faptului ca cele 49 de puncte sunt dispersate în plan (2D) sau în spațiu (3D) iar faptul că ele evoluează independent dar ghidate de efectul de grup va face să se aglomereze chiar în punctul de optim.

Pentru aflarea memoriei ocupate folosim funcția whos din Matlab

Observații generale

În tabelele anterioare au fost sintetizat rezultatele unor indicatori precum durata de funcționare a programului, diferența dintre valoarea optimă și cea rezultată, memoria ocupată, numărul de iterații până la determinarea optimului, etc. În urma acestor rezultate în capitolul următor vor fi trase concluzii și se vor face recomandări pentru folosirea unei metode.

6. Concluzii

Scopul acestei lucrări a fost de a realiza aplicații sugestive pentru metodele de optimizare de tip evolutiv pentru sisteme de mari dimensiuni în vederea comparării acestor metode.

În prima parte a fost definită noțiunea de sistem de mari dimensiuni și problemele care apar:

Număr mare de variabile și de restricții

Complexitatea expresiilor funcției obiectiv și a restricțiilor

Performanțele limitate ale echipamentului de calcul atât din punct de vedere hardware cât și din punct de vedere software

În următoarea etapă au fost prezentate tehnici pentru reducerea unei probleme de mari dimensiuni la subprobleme mai mici dar care sunt interconectare și care sunt caracterizate de existenta restricțiilor. Dintre aceste metode amintim:

tehnici de descompunere

tehnici de partiționare și relaxare

S-a propus un set de metode de tip evolutiv pentru rezolvarea problemelor de optim:

Metode Nelder-Mead (fără restricții)

Metoda BOX (cu restricții)

Algoritmul PSO (Particle Swarm Optimisation)

Pentru cele trei metode de mai sus s-a realizat o prezentare din punct de vedere teoretic, algoritmii au fost implementați în Matlab și s-a realizat un studiu comparativ atât între rezultatele celor trei rutine cât și cu rezultatele obținute prin folosirea unor aplicații realizate de specialiști în domeniu.

Mai exact pentru metodele Nelder – Mead și BOX comparația a fost făcută cu rezultatele obținute pentru aceleași metode dar implementate în produsul software pentru optimizare, SISCON iar pentru PSO comparația s-a făcut cu o aplicație intitulată ToolBox Matlab PSO

Studiul a fost raportat la anumiți indicatori de performanță precum:

timp de rulare

memorie ocupată de algoritm

iterații până la atingerea optimului

eroare

importanța condițiilor inițiale

În funcție de toți acești indicatori se poate observa că fiecare program are puncte slabe și puncte forte. Am notat cu 1 punctajul maxim iau cu 3 punctajul minim.

Pentru exemplele prezentate în lucrare, algoritmul PSO s-a dovedit a fi cel mai bun.

La metoda PSO, deși funcția este evaluată pentru fiecare particulă în parte și în mod normal acest lucru este un inconvenient, faptul ca nu există restricții între ele micșorează volumul de calcul iar rezultatele în ceea ce privește durata de funcționare a programului până la atingerea optimului îl clasează pe primul loc.

Totodată un alt avantaj îl reprezintă posibilitatea de a plasa particulele în mai multe puncte ale domeniului funcției de optimizat iar evoluția lor va fi către mai multe puncte de optim local dar faptul ca întotdeauna se memorează cel mai bun rezultat al grupului, posibilitatea ca optimul obținut să fie și cel global este foarte mare.

Dezavantajul metodelor Nelder-Mead și BOX este că trebuie să rulam programul de mai multe ori cu diferite inițializări și din toate să păstrăm cel mai bun rezultat.

În ceea ce privește numărul de iterații până la atingerea optimului, tot metoda PSO este prima clasată din cauza numărului mare de evaluări ale funcției pentru fiecare iterație.

Algoritmul PSO și-a demonstrat performantele în multe probleme de optimizare. Acest algoritm poate pune probleme în cazul problemelor cu multe puncte de minim local, dar are multiple avantaje față de alte metode.

În implementarea algoritmului PSO nu este nevoie de calcul de gradient, ceea ce face ca implementarea lui să fie mai atractiva fata de algoritmii care cer calculul gradientului. De multe ori este posibil ca funcția noastră să nu poată fi derivabilă și atunci aceste metode evolutive se pretează foarte bine problemei.

Prin modificare parametrilor, se poate modifica viteza de convergenta și performantele algoritmului. De asemenea, impunându-se restricții asupra parametrilor, algoritmul este ușor de controlat.

Având memorie, acest algoritm, în combinație cu algoritmii genetici, poate fi utilizat în procese de învățare.

7. ANEXE – cod algoritmi

7.1. Metoda NELDER – MEAD

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Funcția al cărui minim îl căutam

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function J=optexampfunction(theta)

%Un exemplu de functie de optimizat

J=…

+5*exp(-0.1*((theta(1,1)-15)^2+(theta(2,1)-20)^2))…

-2*exp(-0.08*((theta(1,1)-20)^2+(theta(2,1)-15)^2))…

+3*exp(-0.08*((theta(1,1)-25)^2+(theta(2,1)-10)^2))…

+2*exp(-0.1*((theta(1,1)-10)^2+(theta(2,1)-10)^2))…

-2*exp(-0.5*((theta(1,1)-5)^2+(theta(2,1)-10)^2))…

-4*exp(-0.1*((theta(1,1)-15)^2+(theta(2,1)-5)^2))…

-2*exp(-0.5*((theta(1,1)-8)^2+(theta(2,1)-25)^2))…

-2*exp(-0.5*((theta(1,1)-21)^2+(theta(2,1)-25)^2))…

+2*exp(-0.5*((theta(1,1)-25)^2+(theta(2,1)-16)^2))…

+2*exp(-0.5*((theta(1,1)-5)^2+(theta(2,1)-14

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Nelder-Mead Metoda Simplex (Metoda directa de cautare)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%

%

% Acest program simuleaza minimizarea unei functii folosind

% algoritmul Nelder-Mead

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all %

figure(3)

clf

figure(4)

clf

p=2; % Dimensiunea spatiului de cautare (p+1=numarul de varfuri)

Nds=60; % Numarul maxin de iteratii

beta=1; % Coeficientul de reflexie (ales standard)

gamma=1; % Coeficientul de expansiune(ales standard)

alpha=0.5; % Coeficientul de contractie(ales standard)

% Definirea varfurilor initiale ale simplex-ului:

% Cu urmatoarea initializare, minimul global va fi în punctul (15,5)

%P(:,:,1)=[0 15 30;

% 0 30 0];

% Cu urmatoarea initializare, minimul global va fi în punctul (15,5)

%P(:,:,1)=[0 30 30;

% 15 0 30];

% Cu urmatoarea initializare, minimul global va fi în punctul (15,5)

%P(:,:,1)=[0 0 30;

% 0 30 15];

% Cu urmatoarea initializare, minimul global va fi în punctul (20,15)

% deoarece se pleaca dintr-un punct apropiat de un minim local

%P(:,:,1)=[0 15 30;

% 30 0 30];

% Cu urmatoarea initializare, minimul global va fi în punctul (15,5)

P(:,:,1)=[20 15 15;

5 0 15];

% Cu urmatoarea initializare, minimul global va fi în punctul (20,15)

% deoarece se pleaca dintr-un punct apropiat de un minim local

%P(:,:,1)=[15 20 15;

% 15 15 20];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Realizarea graficului functiei al carui minim în cautam:

x=0:31/100:30; % Pentru functia noastra consideram urmatorul interval în care ia valori

y=x;

% Calcularea functiei pentru care se doreste gasirea minimului

for jj=1:length(x)

for ii=1:length(y)

z(ii,jj)=optexampfunction([x(jj);y(ii)]);

end

end

% Realizam graficul functiei în 3D

figure(1)

clf

surf(x,y,z);

colormap(jet)

% Use next line for generating plots to put în black and white documents.

colormap(gray);

xlabel('x=\theta_1');

ylabel('y=\theta_2');

zlabel('z=J');

title('Functia de minimizat');

% Realizam graficul functiei (curbe de izonivel)

figure(2)

clf

contour(x,y,z,25)

colormap(jet)

% Facem o reprezentare cu tonuri de gri

colormap(gray);

xlabel('x=\theta_1');

ylabel('y=\theta_2');

title('Functia de minimizat (curbe de izonivel)');

hold on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Bucla principala

for j=1:Nds

% Calculam valorile tuturor varfurilor

for i=1:p+1

J(i,j)=optexampfunction([P(1,i,j);P(2,i,j)]);

end

% Cautam cel mai bun și cel mai rau varf și retinem indicii lor

[Jmax(j),maxvertex(j)]=max(J(:,j));

[Jmin(j),minvertex(j)]=min(J(:,j));

thetamax(:,j)=P(:,maxvertex(j),j);

thetamin(:,j)=P(:,minvertex(j),j);

% Cautam centroid-ul tuturor punctelor cu exceptia celui mai rau punct

thetacent(:,j)=0*ones(p,1); % Initializare

for i=1:p+1

thetacent(:,j)=thetacent(:,j)+P(:,i,j);

end

thetacent(:,j)=(1/p)*(thetacent(:,j)-thetamax(:,j));

% Calcularea punctelor de reflectie:

thetaref(:,j)=thetacent(:,j)+beta*(thetacent(:,j)-thetamax(:,j));

Jref(j)=optexampfunction([thetaref(1,j);thetaref(2,j)]);

% Calcularea mai multor valori pentru a realiza testele pentru cel mai bun

% punct reflectat

temp1=J(maxvertex(j),j); % Salvam valoarea varfului care are cel mai rau cost

J(maxvertex(j),j)=-Inf; % Inlocuim elementul cel mai mare cu cel mai mic

[Jmaxwm(j),maxvertexwm(j)]=max(J(:,j)); % Cautam al doilea element ca marime

J(maxvertex(j),j)=temp1; % Aducem matricea la forma initiala

% Urmeaza testele legate de alegerea varfurilor în functie de unele criterii

% de inegalitate

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Primul test la pasul de reflectie

if Jmin(j)>Jref(j),

% Continuam cu expansiunea

thetaexp(:,j)=thetaref(:,j)+gamma*(thetaref(:,j)-thetacent(:,j));

% Urmeaza sa calculam Jexp, care este costul la pasul de expansiune,

Jexp(j)=optexampfunction([thetaexp(1,j);thetaexp(2,j)]);

% Incercam sa realizam expansiunea

if Jexp(j)<Jref(j),

thetanew(:,j)=thetaexp(:,j);

Jnew(j)=Jexp(j);

else

thetanew(:,j)=thetaref(:,j);

Jnew(j)=Jref(j);

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Al doilea test la pasul de reflectie

if Jmaxwm(j)>Jref(j) & Jref(j)>=Jmin(j),

% Punctul de reflectie ia o valoare intermediara

thetanew(:,j)=thetaref(:,j);

Jnew(j)=Jref(j);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Al treilea test la pasul de reflectie

if Jref(j)>=Jmaxwm(j),

% Punctul de reflectie nu e bun

% Realizam retragerea varfului

if Jmax(j)<=Jref(j)

thetanew(:,j)=alpha*thetamax(:,j)+(1-alpha)*thetacent(:,j);

else

thetanew(:,j)=alpha*thetaref(:,j)+(1-alpha)*thetacent(:,j);

end

% Avem nevoie de Jnew pentru a realiza retragerea

Jnew(j)=optexampfunction([thetanew(1,j);thetanew(2,j)]);

end

% Formam noul simplex (utilizam retragerea daca Jnew ne indica faptul ca

% thetanew este mai decat cel mai rau varf:

if Jref(j)>=Jmaxwm(j) & Jnew(j)>Jmax(j),

% Teste legate contractie

% noul varf nu este bun

for i=1:p+1

P(:,i,j+1)=0.5*(P(:,i,j)+thetamin(:,j)); % retragem cel mai bun varf

end

else

% Formam simplexul

% tinem cont ca am folositdeja retragerea

P(:,:,j+1)=P(:,:,j);

P(:,maxvertex(j),j+1)=thetanew(:,j);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% În continuare vor realiza grafice pentru a ilustra cum lucreaza algoritmul

% și ce se intampla la animiti pasi

if j<=4,

figure(3)

subplot(2,2,j)

contour(x,y,z,25)

colormap(jet)

% Am preferat tonuri de gri

colormap(gray);

T=num2str(j-1);

T=strcat('Iteratia j=',T);

ylabel(T)

hold on

% Simplex la iteratia j

simplexplot=plot([P(1,:,j) P(1,:,j)],[P(2,:,j) P(2,:,j)],'bo-');

set(simplexplot,'LineWidth',2);

hold off

%pause

end

%%%%%

if j>4 & j<=8,

figure(4)

subplot(2,2,j-4)

contour(x,y,z,25)

colormap(jet)

% Am preferat tonuri de gri

colormap(gray);

T=num2str(j-1);

T=strcat('Iteratia j=',T);

ylabel(T)

hold on

% Simplex itaratia j

simplexplot=plot([P(1,:,j) P(1,:,j)],[P(2,:,j) P(2,:,j)],'bo-');

set(simplexplot,'LineWidth',2);

hold off

%pause

end

end % Sfarsit bucla principala…

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Cateva grafice legate de rezultatele simularii

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t=0:Nds-1;

figure(5)

clf

contour(x,y,z,25)

colormap(jet)

%Gri

colormap(gray);

xlabel('x=\theta_1');

ylabel('y=\theta_2');

title('Functia de minimizat (curbe de izonivel) ');

hold on

% Adaugam varfurile simplexului care a fost generat

for i=1:p+1,

xx=squeeze(P(1,i,:));

yy=squeeze(P(2,i,:));

vertexplot=plot(xx,yy,'k.');

set(vertexplot,'MarkerSize',10);

end

hold off

% Aprecieri calitative…grafice

figure (6)

clf

subplot(121)

plot(t,thetamin(1,:),'k-',t,thetamin(2,:),'k–')

xlabel('Iteratia j')

ylabel('Pozitia Varfului')

title('Varful al carui cost este minim')

subplot(122)

plot(t,Jmin,'k-',t,Jmax,'k–')

xlabel('Iteratia j')

ylabel('Cost J')

title('Costul celui mai bun și celui mai rau varf')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% End

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

7.2. Metoda BOX

complex_method_main.m

%Algoritmul COMPLEX Program de calcul al punctului de maxim al unei functii neliniare,

%dependente de mai multe variabile. Functia obiectiv trebuie sa aiba variabilele de decizie pozitive.

% Date initiale, ce trebuie definite pentru fiecare problema particulara

clear all

clear figure

global n m k nr_max_iter nr_restr_expl print_control alfa beta gama d iev1 iev2

n=2; %numarul variabilelor explicite

m=2; %numarul restrictiilor explicite

k=6; %numarul punctelor din "complex"

nr_max_iter=500; %numarul maxim de iteratii

nr_restr_expl=2; %numarul restrictiilor implicite

print_control=1; %factor de conditionare a tiparirii:

%print_control=0 tiparirea rezultatului final

%print_control=1 tiparirea rezultatelor intermediare

%plus rezultatul final.

alfa=1.3; %factor de salt

beta=0.0001; %factor de convergenta

gama=5; %iteratii efectuate dupa atingerea convergentei

d=0.001; %factor de corectie

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f=0;

it=0;

iev1=1;

iev2=1;

g=[0 0];

h=[20 20];

c(5)=0;

% Alegerea punctului de start

x(1,1) = 1.0;

x(1,2) = 1.0;

%Generarea numerelor aleatoare V

for u=1:k

for v=1:n

r(u,v) = rand();

end

end

%Deschidem fisierul "complex.ist", ce va con?ine istoria itera?iilor

fp = fopen ('complex.ist', 'w');

fprintf (fp, 'n=%d m=%d k=%d\n', n, m, k);

fprintf (fp, 'nr_max_iter=%d nr_restr__expl=%d alfa=%f \n', nr_max_iter, nr_restr_expl, alfa ) ;

fprintf (fp, 'beta=%f gama=%d delta=%f\n', beta, gama, d );

fprintf (fp, '––––––––––––––––– \n');

for u=1:k

for v=1:n

fprintf (fp, 'r(%d,%d)=%3.4f \n', u, v, r(u,v));

end

end

fprintf (fp, 'numere aleatoare :\n');

[fp, x, r, f, it, iev1, iev2, g, h, c ]=coord (fp, x, r, f, it, iev1, iev2, g, h, c ) ;

if (it > nr_max_iter)

fprintf (fp, 'numarul iteratilor este mai mare de :%3.4f \n', nr_max_iter) ;

exit ( 0);

end

fprintf (fp, '––––––––––––––––– \n');

fprintf (fp,' valoarea finala a functiei :\n');

fprintf (fp,'f[%d]=%3.4f\n', iev2, f(iev2));

fprintf (fp,' valorile finale ale lui x :\n');

for j=1:n

fprintf (fp, 'x[%d][%d]=%3.4f\n', iev2, j, x(iev2,j)) ;

end

fprintf (fp, ' stop !') ;

fclose (fp);

centroid.m

%Subprogramul de calcul al coordonatelor centroidului

function [k1, iev1, c, x, i]=centroid ( k1, iev1, c, x, i)

global n m k nr_max_iter nr_restr_expl print_control alfa beta gama d

rk=0;

for j=1:n

c(j) = 0 ;

for y=1:k

c(j)=c(j) + x(y,j);

end

rk =k1 ;

c(j) = (c(j)-x(iev1,j))/(rk-1);

end

end

functia.m

% Definim functia obiectiv

function [f]=functia( i, x, f)

global n m k nr_max_iter nr_restr_expl print_control alfa beta gama d iev1 iev2

f(i)=-( +5*exp(-0.1*((x(i,1)-15)*(x(i,1)-15)+(x(i,2)-20)*(x(i,2)-20)))…

-2*exp(-0.08*((x(i,1)-20)*(x(i,1)-20)+(x(i,2)-15)*(x(i,2)-15)))…

+3*exp(-0.08*((x(i,1)-25)*(x(i,1)-25)+(x(i,2)-10)*(x(i,2)-10)))…

+2*exp(-0.1 *((x(i,1)-10)*(x(i,1)-10)+(x(i,2)-10)*(x(i,2)-10)))…

-2*exp(-0.5 *((x(i,1)-5 )*(x(i,1)-5 )+(x(i,2)-10)*(x(i,2)-10)))…

-4*exp(-0.1 *((x(i,1)-15)*(x(i,1)-15)+(x(i,2)-5 )*(x(i,2)- 5)))…

-2*exp(-0.5 *((x(i,1)-8 )*(x(i,1)-8 )+(x(i,2)-25)*(x(i,2)-25)))…

-2*exp(-0.5 *((x(i,1)-21)*(x(i,1)-21)+(x(i,2)-25)*(x(i,2)-25)))…

+2*exp(-0.5 *((x(i,1)-25)*(x(i,1)-25)+(x(i,2)-16)*(x(i,2)-16)))…

+2*exp(-0.5 *((x(i,1)-5 )*(x(i,1)- 5)+(x(i,2)-14)*(x(i,2)-14))));

end

restrictii.m

%Subprogramul în care definim valorile limitelor intervalelor de definitie, ale variabilelor

function [g,h,x,i]=restrictii (g, h, x, i)

global n m k nr_max_iter nr_restr_expl print_control alfa beta gama d iev1 iev2

g(1) = 0.0 ;

h(1) = 20.0 ;

g(2) = 0.0 ;

h(2) = 20.0 ;

end

verif.m

%Subprogramul în care verificam daca sunt îndeplinite condi?iile impuse prin restric?ii .

function [ k1, i, g, h, kode, iev1, x, c]=verif ( k1, i, g, h, kode, iev1, x, c)

global n m k nr_max_iter nr_restr_expl print_control alfa beta gama d

e10=1;

while(e10==1)

e10=0;

kt = 0;

[g, h, x, i]=restrictii (g, h, x, i) ;

for j=1:n

if ( x(i,j) < g(j) ) x(i,j) = g(j) + d ;end

if ( x(i,j) > h(j) ) x(i,j) = h(j) – d ;end

if ( x(i,j) < 0 ) x(i,j) = d ;end

end

if ( kode > 0 )

nn = n + 1 ;

for j=nn:m

[g, h, x, i]=restrictii (g, h, x, i) ;

if (x(i,j) < g(j) )

if (x(i,j) > h(j) )

iev1 = i ;

kt = 1 ;

[k1, iev1, c, x, i]=centroid ( k1, iev1, c, x, i);

for l=1:n x(i,l) = ( x(i,l) + c(l)) / 2 ; end

end

end

end

end

if (kt> 0) e10=1; end

end

end

coord.m

%Subrutina de coordonare a tuturor celorlalte subprograme.

%Tipareste toate rezultatele intermediare, în fisierul "complex.ist" .

function [fp, x, r, f, it, iev1, iev2, g, h, c]= coord (fp, x, r, f, it, iev1, iev2, g, h, c)

global n m k nr_max_iter nr_restr_expl print_control alfa beta gama d

i=0; u=0; j=0; k1=0; io=0; kount=0; ia=0; y=0 ;

it = 1 ;

kode=0;

if ( m > n ) kode = 1 ;

end

for u=1:k

for j=1:n

x(u,j) = 0 ;

end

end

fprintf ( fp, 'coordonatele complexului initial:\n') ;

for u=1:k

for j=1:n

i = u ;

[g,h,x,i]=restrictii (g, h, x, i);

x(u,j) = g(j) + r(u,j)*(h(j) – g(j));

end

k1 = u ;

[k1, i, g, h, kode, iev1, x, c]=verif ( k1, i, g, h, kode, iev1, x, c);

if ( print_control == 1 )

if ( u >= 1 )

for j=1:n

fprintf (fp, 'x[%d][%d]=%3.4f\n', u, j, x(u,j));

end

end

end

end

k1 = k;

for i=1:k

[f]=functia( i, x, f);

end

kount = 1 ;

ia = 0 ;%%%%%%

if ( print_control == 1)

fprintf (fp, 'valorile functiei:\n') ;

for j=1:k

fprintf (fp, 'f[%d]=%3.4f\n', j, f(j)) ;

end

end

aux2=1;

while(aux2==1)

aux2=0;

iev1 = 1 ;

for y=2:k

if (f(iev1) > f(y) )

iev1 = y ;

end

end

iev2 = 1 ;

for y=2:k

if (f(iev2) < f(y))

iev2 = y ;

end

end

%Verificarea criteriului de convergenta

aux=0;

if (f(iev2) >= (f(iev1) + beta ))

aux=1;

end

kount = kount + 1 ;

if (( kount < gama )||(aux==1))

%Inocuirea punctului ce genereaz?avaloare mica pentru functie

if (aux==1)

kount=kount-1;

end

if (f(iev2) >= (f(iev1) + beta ))

kount = 1 ;

end

[k1, iev1, c, x, i]=centroid ( k1, iev1, c, x, i) ;

for v=1:n

x(iev1,v) = (1 + alfa ) * c(v) – alfa * x(iev1,v) ;

end

i = iev1 ;

[k1, i, g, h, kode, iev1, x, c]=verif ( k1, i, g, h, kode, iev1, x, c);

[f]=functia( iev1, x, f);

% înlocuirea punctului ce continua sa genereze valoare mica pentru

% functie, cu un punct nou, situat la jumatatea distantei

% dintre punctul vechi și centroid

while(aux2==0)

iev2 = 1 ;

for y=2:k

if (f(iev2) > f(y))

iev2 = y ;

end

end

if (iev2 ~= iev1 )

if ( print_control == 1 )

fprintf (fp,' iteratia nr. %d \n', it) ;

fprintf (fp, ' coord. punctului corectat :\n') ;

fprintf (fp, '––––-iteratia %d ––-

–––––––- \n', it);

for z=1:n

fprintf (fp, 'x[%d][%d]=%3.4f\n',

iev1, z, x(iev1,z)) ;

end

fprintf (fp, '\n');

for i=1:k

fprintf (fp, 'f[%d] = %3.4f\n', i, f(i)) ;

end

fprintf (fp, 'coordonatele centroidului :\n');

for z=1:n

fprintf (fp, 'x_centr[%d] =%3.4f\n', z,c(z)) ;

end

hold on

figure(1)

plot (c(1),c(2),'o');

xlabel('x');

ylabel('y');

title('Evolutia centroidului')

hold off

it = it + 1 ;

if (it <= nr_max_iter)

aux2=1;

end

end

end

if(aux2==0)

for v=1:n

x(iev1,v) = (x(iev1,v) + c(v))/ 2 ;

end

i = iev1 ;

[ k1, i, g, h, kode, iev1, x, c]=

verif ( k1, i, g, h, kode, iev1, x, c ) ;

[f]=functia( i, x, f) ;

end

end

end

end

7.3. PSO (Particle Swarm Optimisation)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Particle Swarm Optimization

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

iterations = 50; %numar de iteratii

inertia = 1; %coeficient de inertie

correction_factor = 2; %factor de corectie

swarm_size = 49; %numarul de puncte

funct_curenta = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Pozitia initiala a grupului

index = 1;

for i = 1 : 7

for j = 1 : 7

swarm(index, 1, 1) = i;

swarm(index, 1, 2) = j;

index = index + 1;

x=i;

y=j;

funct=…

+5*exp(-0.1*((x-15)^2+(y-20)^2))…

-2*exp(-0.08*((x-20)^2+(y-15)^2))…

+3*exp(-0.08*((x-25)^2+(y-10)^2))…

+2*exp(-0.1*((x-10)^2+(y-10)^2))…

-2*exp(-0.5*((x-5)^2+(y-10)^2))…

-4*exp(-0.1*((x-15)^2+(y-5)^2))…

-2*exp(-0.5*((x-8)^2+(y-25)^2))…

-2*exp(-0.5*((x-21)^2+(y-25)^2))…

+2*exp(-0.5*((x-25)^2+(y-16)^2))…

+2*exp(-0.5*((x-5)^2+(y-14)^2));

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = 0;

y = 0;

val = +5*exp(-0.1*((x-15)^2+(y-20)^2))…

-2*exp(-0.08*((x-20)^2+(y-15)^2))…

+3*exp(-0.08*((x-25)^2+(y-10)^2))…

+2*exp(-0.1*((x-10)^2+(y-10)^2))…

-2*exp(-0.5*((x-5)^2+(y-10)^2))…

-4*exp(-0.1*((x-15)^2+(y-5)^2))…

-2*exp(-0.5*((x-8)^2+(y-25)^2))…

-2*exp(-0.5*((x-21)^2+(y-25)^2))…

+2*exp(-0.5*((x-25)^2+(y-16)^2))…

+2*exp(-0.5*((x-5)^2+(y-14)^2));

swarm_best = [min(funct)];

swarm_med = [mean(funct)];

swarm(:, 4, 1) = 1000;% vectorul cu cea mai buna valoare pana acum

swarm(:, 2, 🙂 = 0; % vectorul viteza

gbest=zeros(iterations,1);

for iter = 1 : iterations

funct_curenta = [];

best = swarm_best(end);

%Evaluare pozitie și functia în pozitia curenta

for i = 1 : swarm_size

swarm(i, 1, 1) = swarm(i, 1, 1) + swarm(i, 2, 1)/1.3; %Actualizare pozitia x

swarm(i, 1, 2) = swarm(i, 1, 2) + swarm(i, 2, 2)/1.3; %Actualizare pozitia y

x = swarm(i, 1, 1);

y = swarm(i, 1, 2);

% functia cost

val=…

+5*exp(-0.1*((x-15)^2+(y-20)^2))…

-2*exp(-0.08*((x-20)^2+(y-15)^2))…

+3*exp(-0.08*((x-25)^2+(y-10)^2))…

+2*exp(-0.1*((x-10)^2+(y-10)^2))…

-2*exp(-0.5*((x-5)^2+(y-10)^2))…

-4*exp(-0.1*((x-15)^2+(y-5)^2))…

-2*exp(-0.5*((x-8)^2+(y-25)^2))…

-2*exp(-0.5*((x-21)^2+(y-25)^2))…

+2*exp(-0.5*((x-25)^2+(y-16)^2))…

+2*exp(-0.5*((x-5)^2+(y-14)^2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

funct_curenta = [funct_curenta val];

if val < swarm(i, 4, 1)% daca s-a gasit o pozitie mai buna

swarm(i, 3, 1) = swarm(i, 1, 1); % Actualizare pozitia x

swarm(i, 3, 2) = swarm(i, 1, 2); % Actualizare pozitia y

swarm(i, 4, 1) = val; % Actualizare functie best

end

end

lbest = min(funct_curenta);

swarm_best = [swarm_best min(lbest, best)];

swarm_med = [swarm_med mean(funct_curenta)];

[temp, gbest(iter)] = min(swarm(:, 4, 1)); % Pozitia cea mai buna din grup

%– Actualizare viteza

for i = 1 : swarm_size

swarm(i, 2, 1) = rand*inertia*swarm(i, 2, 1)…

+ correction_factor*rand*(swarm(i, 3, 1)…

– swarm(i, 1, 1))…

+ correction_factor*rand*(swarm(gbest(iter),3,1) …

– swarm(i, 1, 1)); %Componenta x a vitezei

swarm(i, 2, 2) = rand*inertia*swarm(i, 2, 2)…

+ correction_factor*rand*(swarm(i, 3, 2)…

– swarm(i, 1, 2))…

+ correction_factor*rand*(swarm(gbest(iter),3,2)…

– swarm(i, 1, 2)); %Componenta y a vitezei

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Graficul grupului la iteartia curenta

clf

figure(1);

plot(swarm(:, 1, 1), swarm(:, 1, 2), 'x');

title(['Evolutia grupului la iteratia J=',num2str(iter)]);

grid on

axis([-2 30 -2 30]);

pause(.4)

end

figure(2)

stairs([0:1:iterations], swarm_best);

hold on;

stairs([0:1:iterations], swarm_med,'r');

legend('gbest','best mediu');

xlabel ('Nr iteratie');

ylabel ('Valoare functie');

title('Evolutia grupului în cazul inertia = 1 și factorul de corectie = 2.0');

hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%END

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

7.4. ToolBox Matlab PSO

% demopsobehavior.m

% demo of the pso.m function

% the pso tries to find the minimum of the f6 function, a standard

% benchmark

% on the plots, blue is current position, green is Pbest, and red is Gbest

% Brian Birge

% Rev 3.0

% 2/27/08

clear all

close all

clc

help demopsobehavior

warning off

functnames = {'ackley','alpine','DeJong_f2','DeJong_f3','DeJong_f4',…

'Foxhole','Griewank','NDparabola',…

'Rastrigin','Rosenbrock','f6','f6mod','tripod',…

'f6_bubbles_dyn','f6_linear_dyn','f6_spiral_dyn'};

disp('Static test functions, minima don''t change w.r.t. time/iteration:');

disp(' 1) Ackley');

disp(' 2) Alpine');

disp(' 3) DeJong_f2');

disp(' 4) DeJong_f3');

disp(' 5) DeJong_f4');

disp(' 6) Foxhole');

disp(' 7) Griewank');

disp(' 8) NDparabola (for this demo N = 2)');

disp(' 9) Rastrigin');

disp('10) Rosenbrock');

disp('11) Functia mea de test');

disp('12) Schaffer f6 modified (5 f6 functions translated from each other)');

disp('13) Tripod');

disp(' ');

disp('Dynamic test functions, minima/environment evolves over time/iteration:');

disp('14) f6_bubbles_dyn');

disp('15) f6_linear_dyn');

disp('16) f6_spiral_dyn');

functchc=input('Choose test function ? ');

functname = functnames{functchc};

disp(' ');

disp('1) Intense graphics, shows error topology and surfing particles');

disp('2) Default PSO graphing, shows error trend and particle dynamics');

disp('3) no plot, only final output shown, fastest');

plotfcn=input('Choose plotting function ? ');

if plotfcn == 1

plotfcn = 'goplotpso4demo';

shw = 1; % how often to update display

elseif plotfcn == 2

plotfcn = 'goplotpso';

shw = 1; % how often to update display

else

plotfcn = 'goplotpso';

shw = 0; % how often to update display

end

% set flag for 'dynamic function on', only used at very end for tracking plots

dyn_on = 0;

if functchc==15 | functchc == 16 | functchc == 17

dyn_on = 1;

end

%xrng=input('Input search range for X, e.g. [-10,10] ? ');

%yrng=input('Input search range for Y ? ');

xrng=[-30,30];

yrng=[-40,40];

disp(' ');

% if =0 then we look for minimum, =1 then max

disp('0) Minimize')

disp('1) Maximize')

minmax=input('Choose search goal ?');

% minmax=0;

disp(' ');

mvden = input('Max velocity divisor (2 is a good choice) ? ');

disp(' ');

ps = input('How many particles (24 – 30 is common)? ');

disp(' ');

disp('0) Common PSO – with inertia');

disp('1) Trelea model 1');

disp('2) Trelea model 2');

disp('3) Clerc Type 1" – with constriction');

modl = input('Choose PSO model ? ');

% note: if errgoal=NaN then unconstrained min or max is performed

if minmax==1

% errgoal=0.97643183; % max for f6 function (close enough for termination)

errgoal=NaN;

else

% errgoal=0; % min

errgoal=NaN;

end

minx = xrng(1);

maxx = xrng(2);

miny = yrng(1);

maxy = yrng(2);

%–––––––––––––––––––––––––

dims=2;

varrange=[];

mv=[];

for i=1:dims

varrange=[varrange;minx maxx];

mv=[mv;(varrange(i,2)-varrange(i,1))/mvden];

end

ac = [2.1,2.1];% acceleration constants, only used for modl=0

Iwt = [0.9,0.6]; % intertia weights, only used for modl=0

epoch = 400; % max iterations

wt_end = 100; % iterations it takes to go from Iwt(1) to Iwt(2), only for modl=0

errgrad = 1e-99; % lowest error gradient tolerance

errgraditer=100; % max # of epochs without error change >= errgrad

PSOseed = 0; % if=1 then can input particle starting positions, if= 0 then all random

% starting particle positions (first 20 at zero, just for an example)

PSOseedValue = repmat([0],ps-10,1);

psoparams=…

[shw epoch ps ac(1) ac(2) Iwt(1) Iwt(2) …

wt_end errgrad errgraditer errgoal modl PSOseed];

% run pso

% vectorized version

[pso_out,tr,te]=pso_Trelea_vectorized(functname, dims,…

mv, varrange, minmax, psoparams,plotfcn,PSOseedValue);

%–––––––––––––––––––––––––

% display best params, this only makes sense for static functions, for dynamic

% you'd want to see a time history of expected versus optimized global best

% values.

disp(' ');

disp(' ');

disp(['Best fit parameters: ']);

disp([' cost = ',functname,'( [ input1, input2 ] )']);

disp(['–––––––––––']);

disp([' input1 = ',num2str(pso_out(1))]);

disp([' input2 = ',num2str(pso_out(2))]);

disp([' cost = ',num2str(pso_out(3))]);

disp([' mean cost = ',num2str(mean(te))]);

disp([' # of epochs = ',num2str(tr(end))]);

%% optional, save picture

%set(gcf,'InvertHardcopy','off');

%print -dmeta

%print('-djpeg',['demoPSOBehavior.jpg']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%FUNCTIA

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [out]=f6(in)

x=in(:,1);

y=in(:,2);

out=+5.*exp(-0.1.*((x-15).^2+(y-20).^2))…

-2.*exp(-0.08.*((x-20).^2+(y-15).^2))…

+3.*exp(-0.08.*((x-25).^2+(y-10).^2))…

+2.*exp(-0.1.*((x-10).^2+(y-10).^2))…

-2.*exp(-0.5.*((x-5).^2+(y-10).^2))…

-4.*exp(-0.1.*((x-15).^2+(y-5).^2))…

-2.*exp(-0.5.*((x-8).^2+(y-25).^2))…

-2.*exp(-0.5.*((x-21).^2+(y-25).^2))…

+2.*exp(-0.5.*((x-25).^2+(y-16).^2))…

+2.*exp(-0.5.*((x-5).^2+(y-14).^2));

Bibliografie

Dumitru POPESCU, Dan ȘTEFĂNOIU, Ciprian LUPU, Cătălin PETRESCU, Bogdan CIUBOTARU, Călin DIMON Automatică Industrială, Editura Agir, București, 2006

L. LASDON, Optimization Theory for Large Systems, 1970

S. CĂLIN, M. TETIȘCO, I. DUMITRACHE, C.POPEEA, D.POPESCU, Optimizări în Automatizări Industriale, Editura Tehnică, București,1979.

Mircea ANCĂU, Liviu NISTOR Tehnici Numerice de Optimizare în Proiectarea Asistată de Calculator, Editura Tehnica, București, 1996

Vasile NICA, Capitole Speciale ale Cercetărilor Operaționale, Ed. ASE, București, 2001

J. KENNEDY and R. EBERHART, Particle swarm optimization, Dans Proc. of the IEEE Int. Conf. on Neural Networks, Piscataway, NJ, USA, 1995

M. CLERC and J. KENNEDY, The particle swarm-explosion, stability and convergence în a multidimensional complex space. IEEE Trans. Evolutionary Computation, 2002.

J. KENNEDY and W.M. SPEARS, “Matching algorithms to problems: an experimental test of the particle swarm and some genetic algorithms on the multimodal problem generator”. Proc. IEEE Int. Conf. Evolutionary Computation, Anchorage, AK, USA, May 1998.

^2+(y-5).^2))…

-2.*exp(-0.5.*((x-8).^2+(y-25).^2))…

-2.*exp(-0.5.*((x-21).^2+(y-25).^2))…

+2.*exp(-0.5.*((x-25).^2+(y-16).^2))…

+2.*exp(-0.5.*((x-5).^2+(y-14).^2));

Bibliografie

Dumitru POPESCU, Dan ȘTEFĂNOIU, Ciprian LUPU, Cătălin PETRESCU, Bogdan CIUBOTARU, Călin DIMON Automatică Industrială, Editura Agir, București, 2006

L. LASDON, Optimization Theory for Large Systems, 1970

S. CĂLIN, M. TETIȘCO, I. DUMITRACHE, C.POPEEA, D.POPESCU, Optimizări în Automatizări Industriale, Editura Tehnică, București,1979.

Mircea ANCĂU, Liviu NISTOR Tehnici Numerice de Optimizare în Proiectarea Asistată de Calculator, Editura Tehnica, București, 1996

Vasile NICA, Capitole Speciale ale Cercetărilor Operaționale, Ed. ASE, București, 2001

J. KENNEDY and R. EBERHART, Particle swarm optimization, Dans Proc. of the IEEE Int. Conf. on Neural Networks, Piscataway, NJ, USA, 1995

M. CLERC and J. KENNEDY, The particle swarm-explosion, stability and convergence în a multidimensional complex space. IEEE Trans. Evolutionary Computation, 2002.

J. KENNEDY and W.M. SPEARS, “Matching algorithms to problems: an experimental test of the particle swarm and some genetic algorithms on the multimodal problem generator”. Proc. IEEE Int. Conf. Evolutionary Computation, Anchorage, AK, USA, May 1998.

7. ANEXE – cod algoritmi

7.1. Metoda NELDER – MEAD

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Funcția al cărui minim îl căutam

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function J=optexampfunction(theta)

%Un exemplu de functie de optimizat

J=…

+5*exp(-0.1*((theta(1,1)-15)^2+(theta(2,1)-20)^2))…

-2*exp(-0.08*((theta(1,1)-20)^2+(theta(2,1)-15)^2))…

+3*exp(-0.08*((theta(1,1)-25)^2+(theta(2,1)-10)^2))…

+2*exp(-0.1*((theta(1,1)-10)^2+(theta(2,1)-10)^2))…

-2*exp(-0.5*((theta(1,1)-5)^2+(theta(2,1)-10)^2))…

-4*exp(-0.1*((theta(1,1)-15)^2+(theta(2,1)-5)^2))…

-2*exp(-0.5*((theta(1,1)-8)^2+(theta(2,1)-25)^2))…

-2*exp(-0.5*((theta(1,1)-21)^2+(theta(2,1)-25)^2))…

+2*exp(-0.5*((theta(1,1)-25)^2+(theta(2,1)-16)^2))…

+2*exp(-0.5*((theta(1,1)-5)^2+(theta(2,1)-14

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Nelder-Mead Metoda Simplex (Metoda directa de cautare)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%

%

% Acest program simuleaza minimizarea unei functii folosind

% algoritmul Nelder-Mead

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all %

figure(3)

clf

figure(4)

clf

p=2; % Dimensiunea spatiului de cautare (p+1=numarul de varfuri)

Nds=60; % Numarul maxin de iteratii

beta=1; % Coeficientul de reflexie (ales standard)

gamma=1; % Coeficientul de expansiune(ales standard)

alpha=0.5; % Coeficientul de contractie(ales standard)

% Definirea varfurilor initiale ale simplex-ului:

% Cu urmatoarea initializare, minimul global va fi în punctul (15,5)

%P(:,:,1)=[0 15 30;

% 0 30 0];

% Cu urmatoarea initializare, minimul global va fi în punctul (15,5)

%P(:,:,1)=[0 30 30;

% 15 0 30];

% Cu urmatoarea initializare, minimul global va fi în punctul (15,5)

%P(:,:,1)=[0 0 30;

% 0 30 15];

% Cu urmatoarea initializare, minimul global va fi în punctul (20,15)

% deoarece se pleaca dintr-un punct apropiat de un minim local

%P(:,:,1)=[0 15 30;

% 30 0 30];

% Cu urmatoarea initializare, minimul global va fi în punctul (15,5)

P(:,:,1)=[20 15 15;

5 0 15];

% Cu urmatoarea initializare, minimul global va fi în punctul (20,15)

% deoarece se pleaca dintr-un punct apropiat de un minim local

%P(:,:,1)=[15 20 15;

% 15 15 20];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Realizarea graficului functiei al carui minim în cautam:

x=0:31/100:30; % Pentru functia noastra consideram urmatorul interval în care ia valori

y=x;

% Calcularea functiei pentru care se doreste gasirea minimului

for jj=1:length(x)

for ii=1:length(y)

z(ii,jj)=optexampfunction([x(jj);y(ii)]);

end

end

% Realizam graficul functiei în 3D

figure(1)

clf

surf(x,y,z);

colormap(jet)

% Use next line for generating plots to put în black and white documents.

colormap(gray);

xlabel('x=\theta_1');

ylabel('y=\theta_2');

zlabel('z=J');

title('Functia de minimizat');

% Realizam graficul functiei (curbe de izonivel)

figure(2)

clf

contour(x,y,z,25)

colormap(jet)

% Facem o reprezentare cu tonuri de gri

colormap(gray);

xlabel('x=\theta_1');

ylabel('y=\theta_2');

title('Functia de minimizat (curbe de izonivel)');

hold on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Bucla principala

for j=1:Nds

% Calculam valorile tuturor varfurilor

for i=1:p+1

J(i,j)=optexampfunction([P(1,i,j);P(2,i,j)]);

end

% Cautam cel mai bun și cel mai rau varf și retinem indicii lor

[Jmax(j),maxvertex(j)]=max(J(:,j));

[Jmin(j),minvertex(j)]=min(J(:,j));

thetamax(:,j)=P(:,maxvertex(j),j);

thetamin(:,j)=P(:,minvertex(j),j);

% Cautam centroid-ul tuturor punctelor cu exceptia celui mai rau punct

thetacent(:,j)=0*ones(p,1); % Initializare

for i=1:p+1

thetacent(:,j)=thetacent(:,j)+P(:,i,j);

end

thetacent(:,j)=(1/p)*(thetacent(:,j)-thetamax(:,j));

% Calcularea punctelor de reflectie:

thetaref(:,j)=thetacent(:,j)+beta*(thetacent(:,j)-thetamax(:,j));

Jref(j)=optexampfunction([thetaref(1,j);thetaref(2,j)]);

% Calcularea mai multor valori pentru a realiza testele pentru cel mai bun

% punct reflectat

temp1=J(maxvertex(j),j); % Salvam valoarea varfului care are cel mai rau cost

J(maxvertex(j),j)=-Inf; % Inlocuim elementul cel mai mare cu cel mai mic

[Jmaxwm(j),maxvertexwm(j)]=max(J(:,j)); % Cautam al doilea element ca marime

J(maxvertex(j),j)=temp1; % Aducem matricea la forma initiala

% Urmeaza testele legate de alegerea varfurilor în functie de unele criterii

% de inegalitate

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Primul test la pasul de reflectie

if Jmin(j)>Jref(j),

% Continuam cu expansiunea

thetaexp(:,j)=thetaref(:,j)+gamma*(thetaref(:,j)-thetacent(:,j));

% Urmeaza sa calculam Jexp, care este costul la pasul de expansiune,

Jexp(j)=optexampfunction([thetaexp(1,j);thetaexp(2,j)]);

% Incercam sa realizam expansiunea

if Jexp(j)<Jref(j),

thetanew(:,j)=thetaexp(:,j);

Jnew(j)=Jexp(j);

else

thetanew(:,j)=thetaref(:,j);

Jnew(j)=Jref(j);

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Al doilea test la pasul de reflectie

if Jmaxwm(j)>Jref(j) & Jref(j)>=Jmin(j),

% Punctul de reflectie ia o valoare intermediara

thetanew(:,j)=thetaref(:,j);

Jnew(j)=Jref(j);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Al treilea test la pasul de reflectie

if Jref(j)>=Jmaxwm(j),

% Punctul de reflectie nu e bun

% Realizam retragerea varfului

if Jmax(j)<=Jref(j)

thetanew(:,j)=alpha*thetamax(:,j)+(1-alpha)*thetacent(:,j);

else

thetanew(:,j)=alpha*thetaref(:,j)+(1-alpha)*thetacent(:,j);

end

% Avem nevoie de Jnew pentru a realiza retragerea

Jnew(j)=optexampfunction([thetanew(1,j);thetanew(2,j)]);

end

% Formam noul simplex (utilizam retragerea daca Jnew ne indica faptul ca

% thetanew este mai decat cel mai rau varf:

if Jref(j)>=Jmaxwm(j) & Jnew(j)>Jmax(j),

% Teste legate contractie

% noul varf nu este bun

for i=1:p+1

P(:,i,j+1)=0.5*(P(:,i,j)+thetamin(:,j)); % retragem cel mai bun varf

end

else

% Formam simplexul

% tinem cont ca am folositdeja retragerea

P(:,:,j+1)=P(:,:,j);

P(:,maxvertex(j),j+1)=thetanew(:,j);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% În continuare vor realiza grafice pentru a ilustra cum lucreaza algoritmul

% și ce se intampla la animiti pasi

if j<=4,

figure(3)

subplot(2,2,j)

contour(x,y,z,25)

colormap(jet)

% Am preferat tonuri de gri

colormap(gray);

T=num2str(j-1);

T=strcat('Iteratia j=',T);

ylabel(T)

hold on

% Simplex la iteratia j

simplexplot=plot([P(1,:,j) P(1,:,j)],[P(2,:,j) P(2,:,j)],'bo-');

set(simplexplot,'LineWidth',2);

hold off

%pause

end

%%%%%

if j>4 & j<=8,

figure(4)

subplot(2,2,j-4)

contour(x,y,z,25)

colormap(jet)

% Am preferat tonuri de gri

colormap(gray);

T=num2str(j-1);

T=strcat('Iteratia j=',T);

ylabel(T)

hold on

% Simplex itaratia j

simplexplot=plot([P(1,:,j) P(1,:,j)],[P(2,:,j) P(2,:,j)],'bo-');

set(simplexplot,'LineWidth',2);

hold off

%pause

end

end % Sfarsit bucla principala…

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Cateva grafice legate de rezultatele simularii

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t=0:Nds-1;

figure(5)

clf

contour(x,y,z,25)

colormap(jet)

%Gri

colormap(gray);

xlabel('x=\theta_1');

ylabel('y=\theta_2');

title('Functia de minimizat (curbe de izonivel) ');

hold on

% Adaugam varfurile simplexului care a fost generat

for i=1:p+1,

xx=squeeze(P(1,i,:));

yy=squeeze(P(2,i,:));

vertexplot=plot(xx,yy,'k.');

set(vertexplot,'MarkerSize',10);

end

hold off

% Aprecieri calitative…grafice

figure (6)

clf

subplot(121)

plot(t,thetamin(1,:),'k-',t,thetamin(2,:),'k–')

xlabel('Iteratia j')

ylabel('Pozitia Varfului')

title('Varful al carui cost este minim')

subplot(122)

plot(t,Jmin,'k-',t,Jmax,'k–')

xlabel('Iteratia j')

ylabel('Cost J')

title('Costul celui mai bun și celui mai rau varf')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% End

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

7.2. Metoda BOX

complex_method_main.m

%Algoritmul COMPLEX Program de calcul al punctului de maxim al unei functii neliniare,

%dependente de mai multe variabile. Functia obiectiv trebuie sa aiba variabilele de decizie pozitive.

% Date initiale, ce trebuie definite pentru fiecare problema particulara

clear all

clear figure

global n m k nr_max_iter nr_restr_expl print_control alfa beta gama d iev1 iev2

n=2; %numarul variabilelor explicite

m=2; %numarul restrictiilor explicite

k=6; %numarul punctelor din "complex"

nr_max_iter=500; %numarul maxim de iteratii

nr_restr_expl=2; %numarul restrictiilor implicite

print_control=1; %factor de conditionare a tiparirii:

%print_control=0 tiparirea rezultatului final

%print_control=1 tiparirea rezultatelor intermediare

%plus rezultatul final.

alfa=1.3; %factor de salt

beta=0.0001; %factor de convergenta

gama=5; %iteratii efectuate dupa atingerea convergentei

d=0.001; %factor de corectie

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f=0;

it=0;

iev1=1;

iev2=1;

g=[0 0];

h=[20 20];

c(5)=0;

% Alegerea punctului de start

x(1,1) = 1.0;

x(1,2) = 1.0;

%Generarea numerelor aleatoare V

for u=1:k

for v=1:n

r(u,v) = rand();

end

end

%Deschidem fisierul "complex.ist", ce va con?ine istoria itera?iilor

fp = fopen ('complex.ist', 'w');

fprintf (fp, 'n=%d m=%d k=%d\n', n, m, k);

fprintf (fp, 'nr_max_iter=%d nr_restr__expl=%d alfa=%f \n', nr_max_iter, nr_restr_expl, alfa ) ;

fprintf (fp, 'beta=%f gama=%d delta=%f\n', beta, gama, d );

fprintf (fp, '––––––––––––––––– \n');

for u=1:k

for v=1:n

fprintf (fp, 'r(%d,%d)=%3.4f \n', u, v, r(u,v));

end

end

fprintf (fp, 'numere aleatoare :\n');

[fp, x, r, f, it, iev1, iev2, g, h, c ]=coord (fp, x, r, f, it, iev1, iev2, g, h, c ) ;

if (it > nr_max_iter)

fprintf (fp, 'numarul iteratilor este mai mare de :%3.4f \n', nr_max_iter) ;

exit ( 0);

end

fprintf (fp, '––––––––––––––––– \n');

fprintf (fp,' valoarea finala a functiei :\n');

fprintf (fp,'f[%d]=%3.4f\n', iev2, f(iev2));

fprintf (fp,' valorile finale ale lui x :\n');

for j=1:n

fprintf (fp, 'x[%d][%d]=%3.4f\n', iev2, j, x(iev2,j)) ;

end

fprintf (fp, ' stop !') ;

fclose (fp);

centroid.m

%Subprogramul de calcul al coordonatelor centroidului

function [k1, iev1, c, x, i]=centroid ( k1, iev1, c, x, i)

global n m k nr_max_iter nr_restr_expl print_control alfa beta gama d

rk=0;

for j=1:n

c(j) = 0 ;

for y=1:k

c(j)=c(j) + x(y,j);

end

rk =k1 ;

c(j) = (c(j)-x(iev1,j))/(rk-1);

end

end

functia.m

% Definim functia obiectiv

function [f]=functia( i, x, f)

global n m k nr_max_iter nr_restr_expl print_control alfa beta gama d iev1 iev2

f(i)=-( +5*exp(-0.1*((x(i,1)-15)*(x(i,1)-15)+(x(i,2)-20)*(x(i,2)-20)))…

-2*exp(-0.08*((x(i,1)-20)*(x(i,1)-20)+(x(i,2)-15)*(x(i,2)-15)))…

+3*exp(-0.08*((x(i,1)-25)*(x(i,1)-25)+(x(i,2)-10)*(x(i,2)-10)))…

+2*exp(-0.1 *((x(i,1)-10)*(x(i,1)-10)+(x(i,2)-10)*(x(i,2)-10)))…

-2*exp(-0.5 *((x(i,1)-5 )*(x(i,1)-5 )+(x(i,2)-10)*(x(i,2)-10)))…

-4*exp(-0.1 *((x(i,1)-15)*(x(i,1)-15)+(x(i,2)-5 )*(x(i,2)- 5)))…

-2*exp(-0.5 *((x(i,1)-8 )*(x(i,1)-8 )+(x(i,2)-25)*(x(i,2)-25)))…

-2*exp(-0.5 *((x(i,1)-21)*(x(i,1)-21)+(x(i,2)-25)*(x(i,2)-25)))…

+2*exp(-0.5 *((x(i,1)-25)*(x(i,1)-25)+(x(i,2)-16)*(x(i,2)-16)))…

+2*exp(-0.5 *((x(i,1)-5 )*(x(i,1)- 5)+(x(i,2)-14)*(x(i,2)-14))));

end

restrictii.m

%Subprogramul în care definim valorile limitelor intervalelor de definitie, ale variabilelor

function [g,h,x,i]=restrictii (g, h, x, i)

global n m k nr_max_iter nr_restr_expl print_control alfa beta gama d iev1 iev2

g(1) = 0.0 ;

h(1) = 20.0 ;

g(2) = 0.0 ;

h(2) = 20.0 ;

end

verif.m

%Subprogramul în care verificam daca sunt îndeplinite condi?iile impuse prin restric?ii .

function [ k1, i, g, h, kode, iev1, x, c]=verif ( k1, i, g, h, kode, iev1, x, c)

global n m k nr_max_iter nr_restr_expl print_control alfa beta gama d

e10=1;

while(e10==1)

e10=0;

kt = 0;

[g, h, x, i]=restrictii (g, h, x, i) ;

for j=1:n

if ( x(i,j) < g(j) ) x(i,j) = g(j) + d ;end

if ( x(i,j) > h(j) ) x(i,j) = h(j) – d ;end

if ( x(i,j) < 0 ) x(i,j) = d ;end

end

if ( kode > 0 )

nn = n + 1 ;

for j=nn:m

[g, h, x, i]=restrictii (g, h, x, i) ;

if (x(i,j) < g(j) )

if (x(i,j) > h(j) )

iev1 = i ;

kt = 1 ;

[k1, iev1, c, x, i]=centroid ( k1, iev1, c, x, i);

for l=1:n x(i,l) = ( x(i,l) + c(l)) / 2 ; end

end

end

end

end

if (kt> 0) e10=1; end

end

end

coord.m

%Subrutina de coordonare a tuturor celorlalte subprograme.

%Tipareste toate rezultatele intermediare, în fisierul "complex.ist" .

function [fp, x, r, f, it, iev1, iev2, g, h, c]= coord (fp, x, r, f, it, iev1, iev2, g, h, c)

global n m k nr_max_iter nr_restr_expl print_control alfa beta gama d

i=0; u=0; j=0; k1=0; io=0; kount=0; ia=0; y=0 ;

it = 1 ;

kode=0;

if ( m > n ) kode = 1 ;

end

for u=1:k

for j=1:n

x(u,j) = 0 ;

end

end

fprintf ( fp, 'coordonatele complexului initial:\n') ;

for u=1:k

for j=1:n

i = u ;

[g,h,x,i]=restrictii (g, h, x, i);

x(u,j) = g(j) + r(u,j)*(h(j) – g(j));

end

k1 = u ;

[k1, i, g, h, kode, iev1, x, c]=verif ( k1, i, g, h, kode, iev1, x, c);

if ( print_control == 1 )

if ( u >= 1 )

for j=1:n

fprintf (fp, 'x[%d][%d]=%3.4f\n', u, j, x(u,j));

end

end

end

end

k1 = k;

for i=1:k

[f]=functia( i, x, f);

end

kount = 1 ;

ia = 0 ;%%%%%%

if ( print_control == 1)

fprintf (fp, 'valorile functiei:\n') ;

for j=1:k

fprintf (fp, 'f[%d]=%3.4f\n', j, f(j)) ;

end

end

aux2=1;

while(aux2==1)

aux2=0;

iev1 = 1 ;

for y=2:k

if (f(iev1) > f(y) )

iev1 = y ;

end

end

iev2 = 1 ;

for y=2:k

if (f(iev2) < f(y))

iev2 = y ;

end

end

%Verificarea criteriului de convergenta

aux=0;

if (f(iev2) >= (f(iev1) + beta ))

aux=1;

end

kount = kount + 1 ;

if (( kount < gama )||(aux==1))

%Inocuirea punctului ce genereaz?avaloare mica pentru functie

if (aux==1)

kount=kount-1;

end

if (f(iev2) >= (f(iev1) + beta ))

kount = 1 ;

end

[k1, iev1, c, x, i]=centroid ( k1, iev1, c, x, i) ;

for v=1:n

x(iev1,v) = (1 + alfa ) * c(v) – alfa * x(iev1,v) ;

end

i = iev1 ;

[k1, i, g, h, kode, iev1, x, c]=verif ( k1, i, g, h, kode, iev1, x, c);

[f]=functia( iev1, x, f);

% înlocuirea punctului ce continua sa genereze valoare mica pentru

% functie, cu un punct nou, situat la jumatatea distantei

% dintre punctul vechi și centroid

while(aux2==0)

iev2 = 1 ;

for y=2:k

if (f(iev2) > f(y))

iev2 = y ;

end

end

if (iev2 ~= iev1 )

if ( print_control == 1 )

fprintf (fp,' iteratia nr. %d \n', it) ;

fprintf (fp, ' coord. punctului corectat :\n') ;

fprintf (fp, '––––-iteratia %d ––-

–––––––- \n', it);

for z=1:n

fprintf (fp, 'x[%d][%d]=%3.4f\n',

iev1, z, x(iev1,z)) ;

end

fprintf (fp, '\n');

for i=1:k

fprintf (fp, 'f[%d] = %3.4f\n', i, f(i)) ;

end

fprintf (fp, 'coordonatele centroidului :\n');

for z=1:n

fprintf (fp, 'x_centr[%d] =%3.4f\n', z,c(z)) ;

end

hold on

figure(1)

plot (c(1),c(2),'o');

xlabel('x');

ylabel('y');

title('Evolutia centroidului')

hold off

it = it + 1 ;

if (it <= nr_max_iter)

aux2=1;

end

end

end

if(aux2==0)

for v=1:n

x(iev1,v) = (x(iev1,v) + c(v))/ 2 ;

end

i = iev1 ;

[ k1, i, g, h, kode, iev1, x, c]=

verif ( k1, i, g, h, kode, iev1, x, c ) ;

[f]=functia( i, x, f) ;

end

end

end

end

7.3. PSO (Particle Swarm Optimisation)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Particle Swarm Optimization

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

iterations = 50; %numar de iteratii

inertia = 1; %coeficient de inertie

correction_factor = 2; %factor de corectie

swarm_size = 49; %numarul de puncte

funct_curenta = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Pozitia initiala a grupului

index = 1;

for i = 1 : 7

for j = 1 : 7

swarm(index, 1, 1) = i;

swarm(index, 1, 2) = j;

index = index + 1;

x=i;

y=j;

funct=…

+5*exp(-0.1*((x-15)^2+(y-20)^2))…

-2*exp(-0.08*((x-20)^2+(y-15)^2))…

+3*exp(-0.08*((x-25)^2+(y-10)^2))…

+2*exp(-0.1*((x-10)^2+(y-10)^2))…

-2*exp(-0.5*((x-5)^2+(y-10)^2))…

-4*exp(-0.1*((x-15)^2+(y-5)^2))…

-2*exp(-0.5*((x-8)^2+(y-25)^2))…

-2*exp(-0.5*((x-21)^2+(y-25)^2))…

+2*exp(-0.5*((x-25)^2+(y-16)^2))…

+2*exp(-0.5*((x-5)^2+(y-14)^2));

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = 0;

y = 0;

val = +5*exp(-0.1*((x-15)^2+(y-20)^2))…

-2*exp(-0.08*((x-20)^2+(y-15)^2))…

+3*exp(-0.08*((x-25)^2+(y-10)^2))…

+2*exp(-0.1*((x-10)^2+(y-10)^2))…

-2*exp(-0.5*((x-5)^2+(y-10)^2))…

-4*exp(-0.1*((x-15)^2+(y-5)^2))…

-2*exp(-0.5*((x-8)^2+(y-25)^2))…

-2*exp(-0.5*((x-21)^2+(y-25)^2))…

+2*exp(-0.5*((x-25)^2+(y-16)^2))…

+2*exp(-0.5*((x-5)^2+(y-14)^2));

swarm_best = [min(funct)];

swarm_med = [mean(funct)];

swarm(:, 4, 1) = 1000;% vectorul cu cea mai buna valoare pana acum

swarm(:, 2, 🙂 = 0; % vectorul viteza

gbest=zeros(iterations,1);

for iter = 1 : iterations

funct_curenta = [];

best = swarm_best(end);

%Evaluare pozitie și functia în pozitia curenta

for i = 1 : swarm_size

swarm(i, 1, 1) = swarm(i, 1, 1) + swarm(i, 2, 1)/1.3; %Actualizare pozitia x

swarm(i, 1, 2) = swarm(i, 1, 2) + swarm(i, 2, 2)/1.3; %Actualizare pozitia y

x = swarm(i, 1, 1);

y = swarm(i, 1, 2);

% functia cost

val=…

+5*exp(-0.1*((x-15)^2+(y-20)^2))…

-2*exp(-0.08*((x-20)^2+(y-15)^2))…

+3*exp(-0.08*((x-25)^2+(y-10)^2))…

+2*exp(-0.1*((x-10)^2+(y-10)^2))…

-2*exp(-0.5*((x-5)^2+(y-10)^2))…

-4*exp(-0.1*((x-15)^2+(y-5)^2))…

-2*exp(-0.5*((x-8)^2+(y-25)^2))…

-2*exp(-0.5*((x-21)^2+(y-25)^2))…

+2*exp(-0.5*((x-25)^2+(y-16)^2))…

+2*exp(-0.5*((x-5)^2+(y-14)^2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

funct_curenta = [funct_curenta val];

if val < swarm(i, 4, 1)% daca s-a gasit o pozitie mai buna

swarm(i, 3, 1) = swarm(i, 1, 1); % Actualizare pozitia x

swarm(i, 3, 2) = swarm(i, 1, 2); % Actualizare pozitia y

swarm(i, 4, 1) = val; % Actualizare functie best

end

end

lbest = min(funct_curenta);

swarm_best = [swarm_best min(lbest, best)];

swarm_med = [swarm_med mean(funct_curenta)];

[temp, gbest(iter)] = min(swarm(:, 4, 1)); % Pozitia cea mai buna din grup

%– Actualizare viteza

for i = 1 : swarm_size

swarm(i, 2, 1) = rand*inertia*swarm(i, 2, 1)…

+ correction_factor*rand*(swarm(i, 3, 1)…

– swarm(i, 1, 1))…

+ correction_factor*rand*(swarm(gbest(iter),3,1) …

– swarm(i, 1, 1)); %Componenta x a vitezei

swarm(i, 2, 2) = rand*inertia*swarm(i, 2, 2)…

+ correction_factor*rand*(swarm(i, 3, 2)…

– swarm(i, 1, 2))…

+ correction_factor*rand*(swarm(gbest(iter),3,2)…

– swarm(i, 1, 2)); %Componenta y a vitezei

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Graficul grupului la iteartia curenta

clf

figure(1);

plot(swarm(:, 1, 1), swarm(:, 1, 2), 'x');

title(['Evolutia grupului la iteratia J=',num2str(iter)]);

grid on

axis([-2 30 -2 30]);

pause(.4)

end

figure(2)

stairs([0:1:iterations], swarm_best);

hold on;

stairs([0:1:iterations], swarm_med,'r');

legend('gbest','best mediu');

xlabel ('Nr iteratie');

ylabel ('Valoare functie');

title('Evolutia grupului în cazul inertia = 1 și factorul de corectie = 2.0');

hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%END

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

7.4. ToolBox Matlab PSO

% demopsobehavior.m

% demo of the pso.m function

% the pso tries to find the minimum of the f6 function, a standard

% benchmark

% on the plots, blue is current position, green is Pbest, and red is Gbest

% Brian Birge

% Rev 3.0

% 2/27/08

clear all

close all

clc

help demopsobehavior

warning off

functnames = {'ackley','alpine','DeJong_f2','DeJong_f3','DeJong_f4',…

'Foxhole','Griewank','NDparabola',…

'Rastrigin','Rosenbrock','f6','f6mod','tripod',…

'f6_bubbles_dyn','f6_linear_dyn','f6_spiral_dyn'};

disp('Static test functions, minima don''t change w.r.t. time/iteration:');

disp(' 1) Ackley');

disp(' 2) Alpine');

disp(' 3) DeJong_f2');

disp(' 4) DeJong_f3');

disp(' 5) DeJong_f4');

disp(' 6) Foxhole');

disp(' 7) Griewank');

disp(' 8) NDparabola (for this demo N = 2)');

disp(' 9) Rastrigin');

disp('10) Rosenbrock');

disp('11) Functia mea de test');

disp('12) Schaffer f6 modified (5 f6 functions translated from each other)');

disp('13) Tripod');

disp(' ');

disp('Dynamic test functions, minima/environment evolves over time/iteration:');

disp('14) f6_bubbles_dyn');

disp('15) f6_linear_dyn');

disp('16) f6_spiral_dyn');

functchc=input('Choose test function ? ');

functname = functnames{functchc};

disp(' ');

disp('1) Intense graphics, shows error topology and surfing particles');

disp('2) Default PSO graphing, shows error trend and particle dynamics');

disp('3) no plot, only final output shown, fastest');

plotfcn=input('Choose plotting function ? ');

if plotfcn == 1

plotfcn = 'goplotpso4demo';

shw = 1; % how often to update display

elseif plotfcn == 2

plotfcn = 'goplotpso';

shw = 1; % how often to update display

else

plotfcn = 'goplotpso';

shw = 0; % how often to update display

end

% set flag for 'dynamic function on', only used at very end for tracking plots

dyn_on = 0;

if functchc==15 | functchc == 16 | functchc == 17

dyn_on = 1;

end

%xrng=input('Input search range for X, e.g. [-10,10] ? ');

%yrng=input('Input search range for Y ? ');

xrng=[-30,30];

yrng=[-40,40];

disp(' ');

% if =0 then we look for minimum, =1 then max

disp('0) Minimize')

disp('1) Maximize')

minmax=input('Choose search goal ?');

% minmax=0;

disp(' ');

mvden = input('Max velocity divisor (2 is a good choice) ? ');

disp(' ');

ps = input('How many particles (24 – 30 is common)? ');

disp(' ');

disp('0) Common PSO – with inertia');

disp('1) Trelea model 1');

disp('2) Trelea model 2');

disp('3) Clerc Type 1" – with constriction');

modl = input('Choose PSO model ? ');

% note: if errgoal=NaN then unconstrained min or max is performed

if minmax==1

% errgoal=0.97643183; % max for f6 function (close enough for termination)

errgoal=NaN;

else

% errgoal=0; % min

errgoal=NaN;

end

minx = xrng(1);

maxx = xrng(2);

miny = yrng(1);

maxy = yrng(2);

%–––––––––––––––––––––––––

dims=2;

varrange=[];

mv=[];

for i=1:dims

varrange=[varrange;minx maxx];

mv=[mv;(varrange(i,2)-varrange(i,1))/mvden];

end

ac = [2.1,2.1];% acceleration constants, only used for modl=0

Iwt = [0.9,0.6]; % intertia weights, only used for modl=0

epoch = 400; % max iterations

wt_end = 100; % iterations it takes to go from Iwt(1) to Iwt(2), only for modl=0

errgrad = 1e-99; % lowest error gradient tolerance

errgraditer=100; % max # of epochs without error change >= errgrad

PSOseed = 0; % if=1 then can input particle starting positions, if= 0 then all random

% starting particle positions (first 20 at zero, just for an example)

PSOseedValue = repmat([0],ps-10,1);

psoparams=…

[shw epoch ps ac(1) ac(2) Iwt(1) Iwt(2) …

wt_end errgrad errgraditer errgoal modl PSOseed];

% run pso

% vectorized version

[pso_out,tr,te]=pso_Trelea_vectorized(functname, dims,…

mv, varrange, minmax, psoparams,plotfcn,PSOseedValue);

%–––––––––––––––––––––––––

% display best params, this only makes sense for static functions, for dynamic

% you'd want to see a time history of expected versus optimized global best

% values.

disp(' ');

disp(' ');

disp(['Best fit parameters: ']);

disp([' cost = ',functname,'( [ input1, input2 ] )']);

disp(['–––––––––––']);

disp([' input1 = ',num2str(pso_out(1))]);

disp([' input2 = ',num2str(pso_out(2))]);

disp([' cost = ',num2str(pso_out(3))]);

disp([' mean cost = ',num2str(mean(te))]);

disp([' # of epochs = ',num2str(tr(end))]);

%% optional, save picture

%set(gcf,'InvertHardcopy','off');

%print -dmeta

%print('-djpeg',['demoPSOBehavior.jpg']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%FUNCTIA

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [out]=f6(in)

x=in(:,1);

y=in(:,2);

out=+5.*exp(-0.1.*((x-15).^2+(y-20).^2))…

-2.*exp(-0.08.*((x-20).^2+(y-15).^2))…

+3.*exp(-0.08.*((x-25).^2+(y-10).^2))…

+2.*exp(-0.1.*((x-10).^2+(y-10).^2))…

-2.*exp(-0.5.*((x-5).^2+(y-10).^2))…

-4.*exp(-0.1.*((x-15).^2+(y-5).^2))…

-2.*exp(-0.5.*((x-8).^2+(y-25).^2))…

-2.*exp(-0.5.*((x-21).^2+(y-25).^2))…

+2.*exp(-0.5.*((x-25).^2+(y-16).^2))…

+2.*exp(-0.5.*((x-5).^2+(y-14).^2));