7. Rozwiązywanie zadań opisanych równaniami nieliniowymi Pierwiastki rzeczywiste równania nieliniowego o ogólnej postaci f(x)=0 można wyznaczyć za pomocą funkcji fzero.
Tekst
funkcji
fzero można wyświetlić na ekranie – podobnie jak wiele innych funkcji - pisząc polecenie
type
fzero
Opis wywołania funkcji fzero jest dostępny za pomocą polecenia help
fzero
Wywołanie funkcji z 2 argumentami wejściowymi
p = fzero(‘nazwafunkcji’,x0)
daje w wyniku pierwiastek równania znajdujący się w otoczeniu punktu x0.
Przykładowo
polecenie
p = fzero(‘sin’,3)
daje w wyniku wartość pierwiastka p = pi, co odpowiada przejściu funkcji sinus przez zero w punkcie 180O .
Wyliczone pierwiastki równania nieliniowego f(x)=0 można pokazać na wykresie za pomocą funkcji fplot.
Jeśli wykładniki potęg zmiennej x w równaniu nieliniowym f(x)=0 są liczbami naturalnymi, to do rozwiązania zaleca się stosowanie funkcji roots wyznaczającej miejsca zerowe wielomianu.
Polecenie
p = roots(c)
pozwala wyznaczyć pierwiastki wielomianu o następującej postaci n
n 1
−
2
(
y x ) = c x + c x
+ ... + c x + c x + c (7.1) 1
2
n 1
−
n
n 1
+
Należy zwrócić uwagę, że postać wielomianu podana wzorem (7.1) różni się od postaci stosowanej w podręcznikach matematyki, w których zwykle wyraz wolny ma indeks 0.
73
Przykład
rozwiązania równania nieliniowego za pomocą fzero oraz roots Należy wyznaczyć wszystkie rozwiązania równania nieliniowego x3 + x2 − 3x − 3 = 0 (7.2) za pomocą funkcji fzero oraz roots, a następnie zilustrować na wykresie funkcji położenie tych rozwiązań.
function y=wielom(x)
% wielomian, w ktorym zmienna x moze byc wektorem
y=x.^3 + x.^2 -3.*x -3; return
function yzero=osx(x)
% zerowa os odcietych
yzero=0; return
function [xroots,xfzero] = oblicz
% funkcja znajduje pierwiastki za pomoca roots oraz fzero
% wielomian: y = x^3+x^2-3*x-3
format long; % zwiekszenie ilosci wyswietlanych cyfr
disp('pierwiastki wyznaczone za pomoca roots :');
xroots = roots([1 1 -3 -3]);
disp(xroots);
disp(' sprawdzenie dokladnosci za pomoca funkcji polyval'); epsxroots = polyval([1 1 -3 -3], xroots); disp(epsxroots); disp('pierwiastki wyznaczone za pomoca fzero :');
disp('pierwiastek 1-szy w otoczeniu x0=2 :');
x1=fzero('wielom',2); disp(x1);
disp('pierwiastek 2-gi w otoczeniu x0=-2 :');
x2=fzero('wielom',-2); disp(x2);
disp('pierwiastek 3-ci w otoczeniu x0=0 :');
x3=fzero('wielom',0); disp(x3);
xfzero=[x1; x2; x3];
disp('roznica rozwiazan xroots oraz xfzero');
roznica=xroots-xfzero; disp(roznica);
% Wykres wielomianu y=x^3+x^2-3*x-3
fplot('[wielom(x),osx(x)]',[-2,2]);
hold on;
plot(x1,0,'kO',x2,0,'bO',x3,0,'rO');
text(x1,0.5,'x1'); text(x2,0.5,'x2'); text(x3,0.5,'x3'); title(' Wielomian y = x^3 +x^2 - 3*x -3 ');
xlabel(' x'); ylabel(' y'); grid on;
return
Po wywołaniu funkcji oblicz otrzymuje się wykres pokazany na rys. 7.2. oraz wyniki obliczeń w oknie poleceń Matlaba.
Wielomian y = x3 + x2 - 3x - 3
3
2
1
x2
y
x3
x1
0
-1
-2
-3
-4
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
x
Rys. 7.1. Wykres wielomianu z naniesionymi wartościami pierwiastków
>> [w,z]=oblicz
pierwiastki wyznaczone za pomoca roots :
1.73205080756888
-1.73205080756888
-1.00000000000000
sprawdzenie dokladnosci za pomoca funkcji polyval
1.0e-014 *
0.22204460492503
-0.39968028886506
0
pierwiastki wyznaczone za pomoca fzero :
pierwiastek 1-szy w otoczeniu x0=2 :
1.73205080756888
pierwiastek 2-gi w otoczeniu x0=-2 :
-1.73205080756888
pierwiastek 3-ci w otoczeniu x0=0 :
75
-1.00000000000000
roznica rozwiazan xroots oraz xfzero
1.0e-014 *
0
-0.11102230246252
0
w =
1.73205080756888
-1.73205080756888
-1.00000000000000
z =
1.73205080756888
-1.73205080756888
W nowszych wersjach Matlaba funkcję fzero można stosować do rozwiązywania równania nieliniowego zdefiniowanego w postaci funkcji z parametrami wejściowymi.
Zadania elektrotechniczne są często opisane układami równań nieliniowych. Do ich rozwiązania można wykorzystać funkcję fsolve.
Układ
równań nieliniowych ma następującą ogólną postać
f(x) = 0 (7.3)
gdzie
f(x) = [f (x), f (x), ... f (x)] (7.4)
1
2
n
x = [ x ,x ,...,x ] (7.5) 1
2
n
Funkcja
fsolve rozwiązuje układ równań nieliniowych metodą najmniejszych kwadratów. Należy zwrócić uwagę, że funkcja fsolve w nowszych wersjach Matlaba korzysta z tablic strukturalnych, w których dostęp do danych odbywa się przez podanie nazwy pola. Tablice strukturalne mogą zawierać różne typy danych, np. jedno pole zawiera liczbę, inne łańcuch.
Tablicę strukturalną można utworzyć za pomocą instrukcji przypisania albo za pomocą funkcji struct:
• za pomocą instrukcji przypisania, np.
tabstrukt.nazwa=’pomiar 1’;
tabstrukt.srednia=37.31;
tabstrukt.pomiar=[1 2 3 4 5 6 7 8 9];
• za pomocą funkcji tabstruct=struct(‘pole1’,wartosc1,’pole2’,wartosc2,..) tabstrukt=...
struct('nazwa','pomiar 1','srednia',57,'pomiar',[1 2 3 4 5]) tabstrukt =
nazwa: 'pomiar 1'
srednia: 57
pomiar: [1 2 3 4 5]
Aby
powiększyć tablicę strukturalną, dodaje się indeks za nazwą tablicy przed polem
tabstrukt(2).nazwa=’pomiar 2’;
tabstrukt(2).srednia=17.31;
tabstrukt(2).pomiar=[11 12 13 14 15 16 17 18 19];
Podczas rozszerzania struktury Matlab wypełnia niewyspecyfikowane pola macierzami pustymi. Dostęp do danego pola jest możliwy po podaniu nazwy struktury i nazwy pola
tabstrukt(1).pomiar
ans =
1 2 3 4 5
Dostęp do danego elementu wewnątrz pola jest możliwy po przez podanie indeksów
tabstrukt(1).pomiar(1,3)
ans =
3
Inny sposób dostępu wykorzystuje funkcję getfield
f=getfield(nazwa_tablicy,{ineks_tablicy},’nazwa_pola’,{indeks_pola} ) Pola
są usuwane za pomocą funkcji rmfiled
rmfiled(nazwa_tablicy,’nazwa_pola’)
Nazwy pól tablicy strukturalnej można odczytać za pomocą funkcji fieldnames.
Przykładowo
tabstrukt.nazwa='pomiar 1';
tabstrukt.srednia=3.7;
tabstrukt.pomiar=[1 2 3 4 5 6 7 8 9];
fieldnames(tabstrukt)
ans =
77
'nazwa'
'srednia'
'pomiar'
Warianty
wywołań funkcji fsolve
1. x = fsolve(fun,x0)
gdzie fun – funkcja akceptująca parametr wejściowy x, zwracająca wartość f(x).
2. x = fsolve(fun,x0,options)
gdzie options – opcje poszukiwania rozwiązania ustalane poprzez wywołanie funkcji foptions lub w nowszych wersjach Matlaba przez optimset.
3. x = fsolve(fun,x0,[])
gdzie [] – pusta macierz zamiast wektora options powoduje przyjęcie standardowych wartości opcji.
4. x = fsolve(fun,x0,options,p1,p2,...)
lub
x = fsolve(fun,x0,[],p1,p2,...)
gdzie p1,p2,... – parametry funkcji fun od których zależy rozwiązanie.
5. [x,fval] = fsolve(fun,x0,options,p1,p2,...)
gdzie fval – wartość funkcji celu, czyli sumy kwadratów błędów rozwiązań.
6. [x,fval,exitflag] = fsolve(fun,x0,options,p1,p2,...) gdzie
exitflag >0 – rozwiązanie zbieżne,
exitflag = 0 – osiągnięto maksymalną liczbę iteracji, exitflag < 0 – brak zbieżności procesu iteracyjnego.
7. [x,fval,exitflag,output] = fsolve(fun,x0,options,p1,p2,...) gdzie
output.iterations – liczba iteracji,
output.funcCount – liczba szacowania funkcji celu, output.algorithm – wykorzystany algorithm,
output.cgiterations – liczba iteracji CG,
output.firstorderopt – optymalizacja rzędu pierwszego.
8. [x,fval,exitflag,output,jacob] = fsolve(fun,x0,options,p1,p2,...) gdzie jacob – macierz Jacobiego funkcji f(x) w punkcie x.
W Matlabie 4.2 funkcję fsolve wywołuje się z podaniem gradientu stosując następującą postać polecenia
options=foptions;
x=fsolve(fun,x0,options,grad,p1,p2,..);
gdzie
grad – m-plik zawierający gradient funkcji y=f(x).
Należy zwrócić uwagę na fakt, że macierze w Matlabie zapamiętywane są kolejnymi kolumnami, a nie wierszami jakby się wydawało. Z tego powodu gradient jest równy transponowanej macierzy Jacobiego danego układu równań nieliniowych.
Wykorzystanie
funkcji
fsolve do rozwiązywania układu równań nieliniowych zilustrowano na przykładzie równań węzłowych napięciowo-mocowych.
Przykład zastosowania funkcji fsolve do analizy układu przesyłowego: system-linia-odbiór
Duży zakład przemysłowy jest zasilany linią wysokiego napięcia 220 kV, rys. 7.2.
Moc czynna i bierna dopływająca do odbioru wyrażona jest układem 2 równań węzłowych kwadratowych, zapisanych w jednostkach względnych odniesionych do mocy bazowej 484 MV⋅A oraz napięcia znamionowego 220 kV
− P = GU 2 − Ge − Bf (7.6)
− Q = − BU 2 + Be − Gf (7.7) gdzie
P, Q – moc czynna i bierna odbioru,
G, B – konduktancja i susceptancja podłużna linii, U = e + jf – napięcie zespolone odbioru zapisane w postaci algebraicznej, e – składowa prostokątna rzeczywista napięcia odbioru, f – składowa prostokątna urojona napięcia odbioru, U2 = e2 + f2 - kwadrat modułu napięcia odbioru, U =
2
2
e + f
– wartość skuteczna napięcia odbioru.
Należy wyznaczyć napięcie odbioru odpowiadające poborowi mocy czynnej i biernej. Dodatkowo należy wyznaczyć zależność napięcia odbioru od zmieniającego się poboru mocy biernej, przy stałym poborze mocy czynnej.
79
Z=R+jX S=P+jQ
SEE
U
s
U
Im
Us=Us+j0
e
Re
U=e+jf
f
Rys. 7.2. Schemat układu zasilającego duży zakład przemysłowy oraz wykres napięć w układzie.
Równania
węzłowe opisujące układ przesyłowy pokazany na rys. 7.2 mają postać Ge2 + Gf 2 − Ge − Bf + P = 0 (7.8)
− Be2 − Bf 2 + Be − Gf + Q = 0 (7.9) Układ
równań napięciowo-mocowych można rozwiązać za pomocą fsolve. Przy pisaniu programu wygodnie jest dokonać podstawień pozwalających wykorzystać standardowe sposoby wywołania funkcji fsolve .
W naszym przypadku będą to następujące podstawienia
x(1)=e, x(2)=f, p1= G, p2= B, p3=P, p4=Q
W celu rozwiązania zadania za pomocą funkcji fsolve utworzono 3 następujące funkcje:
• rmf() – funkcja definiująca układ równań nieliniowych,
• rmg() – funkcja obliczająca gradient równań węzłowych,
• rmobl() – funkcja sterująca obliczeniami.
function f = rmf(x,G,B,P,Q)
% rownania wezlowe mocy
f(1)= G*x(1)^2 + G*x(2)^2 - G*x(1) - B*x(2) + P;
f(2)=-B*x(1)^2 - B*x(2)^2 + B*x(1) - G*x(2) + Q;
return
function df = rmg(x,G,B,P,Q)
% gradient rownan wezlowych mocy odbioru
% f(1)= G*x(1)^2 + G*x(2)^2 - G*x(1) - B*x(2) + P;
% f(2)=-B*x(1)^2 - B*x(2)^2 + B*x(1) - G*x(2) + Q;
macJac=[ 2*G*x(1)-G
2*G*x(2)-B
-2*B*x(1)+B
-2*B*x(2)-G];
df = macJac'; %gradient rowny transponowanej m. Jacobiego return
function [Q,U] = rmobl(R,X,Podb,Qodb,Sb,Un)
% function [Q,U] = rmroz(P,Q,Sb,Un)
% wyznacza charakterystyke Q-U mocy biernej
% dostarczanej linia wysokiego napiecia
% R, X - rezystancja i reaktancja linii w omach,
% P - moc czynna odbioru, MW
% Q - moc bierna odbioru, Mvar
% Sb - moc bazowa, MVA
% Un - napiecie znamionowe linii, kV
% Q - wektor mocy biernej zmienianej od 0 do Qmax, Mvar
% U - wartosc skuteczna napiecia odbioru, kV
if nargin < 6 Un=220; end
if nargin < 5 Sb=484; end
if nargin < 4 Podb=242; end
if nargin < 3 Qodb=121; end
if nargin < 2 R=10; end
if nargin < 1 X=20; end
j=sqrt(-1);
Zb=Un^2/Sb; % impedancja bazowa
Z=(R+j*X)/Zb; Y=1/Z;
G=real(Y); B=imag(Y); % parametry linii w jedn. wzglednych P=Podb/Sb; Q=Qodb/Sb; % przeliczenie mocy na jednostki wzgledne x0=[1 0]; x=x0; % punkt startowy obliczen
options=foptions; % opcje standardowe
% Wyznaczenie napiecia odbioru
x=fsolve('rmf',x0,options,'rmg',G,B,P,Q);
u=sqrt(x(1).^2 + x(2)^2); % wyznaczone napiecie odbioru Urm = u*Un; % napiecie odbioru w kV
Qrm=Q*Sb; % moc bierna odbioru w Mvar
disp('Napiecie odbioru w kV wynosi:');
disp(Urm);
% Wyznaczanie charakterystyki napieciowej Q-U
i=0; Q=0;
U=[]; u0=1; U(1)=u0; u=1;
Qodb=[]; Qodb(1)=Q;
Ulawina=0.65*u0; % przyblizony punkt lawiny napiecia
81
while u > Ulawina % obliczenia az do lawiny napiecia i=i+1; % kolejny krok charakterystyki Q-U
Q=Q+0.01; % przyrost mocy biernej
x=fsolve('rmf',x0,options,'rmg',G,B,P,Q);
u=sqrt(x(1).^2 + x(2)^2); % napiecie
Qodb(i)=Q*Sb; % moc bierna w Mvar
Uodb(i)=u*Un; % napiecie w kV
Urmi(i)=Urm; % napiecie odbioru
end
nroz=i; Qlawina=Q*Sb; Ulawina=Ulawina*Un;
disp(' Liczba wyznaczonych punktow charakterystyki
napieciowej');
disp(nroz);
plot(Qodb,Uodb,'b-',Qodb,Urmi,'r--');
title('Charakterystyka Q-U mocy biernej dostarczanej linia 220
kV');
xlabel('Q, Mvar'); ylabel('U, kV');
legend('Krzywa Q-U','Napiecie odbioru');
hold on;
plot(Qrm,Urm,'O');
text(Qrm,Urm+2,'Punkt pracy ukladu przesylowego');
plot(Qlawina,Ulawina,'O');
text(Qlawina-150,Ulawina+1,'Punkt lawiny napiecia');
grid on;
disp('KONIEC');
return
Przebieg
obliczeń numerycznych po wywołaniu funkcji rmobl jest następujący Napiecie odbioru w kV wynosi:
194.2949
Liczba wyznaczonych punktow charakterystyki napieciowej 89
KONIEC
Funkcja
rmobl() po wyznaczeniu kolejnych rozwiązań za pomocą funkcji fsolve tworzy wykres zmian napięcia odbioru od poboru mocy biernej pokazany na rys. 7.3.
Charakterystyka Q-U mocy biernej dostarczanej linia 220 kV
210
Krzywa Q-U
Napiecie odbioru
200
Punkt pracy ukladu przesylowego
190
180
U,
kV
170
160
150
Punkt lawiny napiecia
140 0
50
100
150
200
250
300
350
400
450
Q, Mvar
Rys. 7.3. Wyniki analizy układu przesyłowego: system-linia-odbiór otrzymane za pomocą funkcji solve
Zadania do samodzielnego rozwiązania
Zadanie
1
Zmodyfikować funkcje rmf(), rmg(), rmobl() tak, aby uzyskać charakterystykę zmian napięcia od zmian przesyłanej mocy czynnej zamiast zmian mocy biernej.
Zadanie
2
Wzorując się na funkcjach rmf(), rmg(), rmobl() opracować program rozwiązujący następujący układ równań nieliniowych
x x + bx − bx − b2 = 0 (7.10) 1 2
1
2
x x + bx − 2bx − 2b2 = 0 (7.11) 1 2
1
2
Parametr b może przyjmować dowolną wartość z przedziału 0.2 ≤ b ≤ 0.8 (7.12)
Przyjąć punkt startowy x0=[2*b+0.1 b+0.1];