WYDZIAŁ ETI PG
Katedra Systemó w Elektróniki Mórskiej
Technika obliczeniowa i symulacyjna
- laboratorium
Automatyczne konstruowanie równań liniowych
opisujących obwód elektryczny
oraz ich numeryczne rozwiązywanie.
MATLAB i SPICE jako narzędzia do obliczania prądów, napięć
i mocy oraz elementów zastępczych w obwodzie elektrycznym
(Materiały dó ćwiczenia nr 6)
Opracówał: Czesław Stefan ski
Gdan sk 2015
C. Stefański, Materiały do ćw 6 Lab. TOiS
2/12
Ostatnia modyfikacja: maj 2016
1. Wstęp
Niniejsze ćwiczenie jest kolejnym, które jest poświęcone wdrożeniu umiejętności obliczania prądów i napięć w ob-
wodzie przy wykorzystaniu narzędzia MATLAB z jego algorytmami numerycznego rozwiązywania układów równań linio-
wych oraz nabyciu obycia w korzystaniu z symulacyjnego narzędzia Micro̵CAP (jedna z „odmian” SPICE’a)
1
. Jednak tym
razem będziemy starali się tak zorganizować proces badania obwodu, by jeszcze bardziej odciążyć użytkownika od „ręcz-
nego” udziału w tym badaniu – w tym celu spróbujemy „zmusić” komputer, by nie tylko rozwiązywał układy równań
opisujące obwód, ale by je również sam tworzył! W tej sytuacji warto przypomnieć, że przy tworzeniu algorytmów
analizy automatycznej (komputerowej) wygodnie jest trzymać się zasady, że każdy dwójnik ma „swój” prąd i „swoje”
napięcie, przeciwnie względem siebie zastrzałkowane (czyli stosujemy tzw. strzałkowanie skojarzone). Przez to algo-
rytmy stają się bardziej jednorodne i dzięki temu zwykle są szybsze mimo większej niż przy strzałkowaniu „chaotycznym”
liczby niewiadomych.
Przypominamy, że na rysunkach 1 mamy do czynienia z tzw. skojarzonym strzałkowaniem prądów i napięć – polega
ono na tym, że zwrot spadku napięcia na każdym z elementów obwodu jest przeciwny do kierunku przepływającego
przez ten element prądu.
W prostym obwodzie elektrycznym, na przykład składającym się z podobnej liczby elementów jak w obwodzie
przedstawionym na rys. 1a, poszukiwanie takich podstawowych wielkości elektrycznych jak natężenie prądu (prąd) 𝒊,
czy też różnica potencjałów (napięcie) 𝒖, może być względnie łatwym zadaniem, również dlatego, że są to liniowe ob-
wody prądu stałego. Jednak ze wzrostem liczby i różnorodności elementów tworzących strukturę obwodu stopniowo
zmniejsza się możliwość – czy też opłacalność – ręcznego (analitycznego) poszukiwania prądów i napięć gałęziowych.
i
0
i
1
i
2
i
3
E
0
R
1
R
2
R
3
u
0
u
2
u
3
u
1
J
5
u
5
i
5
i
1
i
2
i
3
i
4
E
1
R
2
R
3
R
4
u
1
u
3
u
4
u
2
t
Rys. 1a. Prosty obwód prądu stałego
Rys. 1b. Przykład strzałkowania skojarzonego
2. Poszukiwanie prądów, napięć i mocy w rozbudowanym obwodzie elektrycznym
Dla ustalenia uwagi przyjmiemy, że naszym zadaniem będzie dokonanie w miarę automatycznej i wszechstronnej
analizy obwodu z rysunku 1b, w tym wyznaczenie prądów i napięć, a także mocy na każdym z jego elementów. Ponadto
porównamy wyniki naszej analizy z wynikami symulacji przeprowadzonej w Micro
̵
Cap-ie, a także sprawdzimy bilanso-
wanie się mocy.
Jednak tym razem wykorzystamy opis obwodu elementarnego poziomu, tzw. opis tableau, którego tworzenie jest
w miarę łatwe do oprogramowania.
Pozostałe pojęcia, algorytmy i narzędzia, wykorzystywane w tym ćwiczeniu, zostały już opisane w materiale do
wcześniejszych ćwiczeń. W szczególności wykorzystamy macierz incydencji 𝐴 grafu obwodu oraz tzw. transformację
węzłową 𝑼 = 𝑨
𝑻
⋅ 𝑽.
.
.
.
.
6
3
7
4
1
2
5
9
8
k
l
i
"
A
= [
1 -1
1
0 0
0
0 0
0
0
1
0 -1 1
0 -1 0
0
0
0
0
0 0 -1
1 1 -1
0
0 -1
1 0
1
0 0
1
]
Rys. 2. Przykładowy graf prądowy i jego zredukowana macierz incydencji
1
Narzędzie to, po „narysowaniu” obwodu w jego oknie graficznym przez użytkownika, potrafi automatycznie, acz bez ujawniania ich użytkownikowi,
utworzyć odpowiednie równania opisujące obwód, rozwiązać je (podobnymi jak MATLAB algorytmami) i udostępnić rozwiązania międ zy innymi we
wspomnianym oknie z obwodem.
C. Stefański, Materiały do ćw 6 Lab. TOiS
3/12
Ostatnia modyfikacja: maj 2016
Graf z rysunku 2, przy założeniu, że węzeł i obwodu jest węzłem odniesienia, stanowił podstawę do zapisania jego
zredukowanej macierzy incydencji 𝑨 = [𝒂
𝒌𝒍
] (dalej nazywanej krótko macierzą incydencji) według następującej umowy:
𝑎
𝑘𝑙
= 1 w przypadku, gdy prąd w gałęzi o numerze „𝑙” wypływa z węzła „𝑘”,
𝑎
𝑘𝑙
= -1 w przypadku, gdy prąd w gałęzi o numerze „𝑙” wpływa do węzła „𝑘”,
𝑎
𝑘𝑙
= 0 w przypadku, gdy gałąź o numerze „𝑙” nie łączy się z węzłem „𝑘”.
Metoda tableau opisu obwodu
Dla obwodu o 𝑏 gałęziach oraz (𝑛 + 1) węzłach [dla obwodu z rysunku 1b mamy 𝑏 = 5 i 𝑛 = 2] węzły numerujemy
tak, by węzeł odniesienia (masy) miał numer zero, a pozostałe węzły numery od 1 do 𝑛. Wtedy napięcia i prądy w takich
gałęziach można przedstawić w formie następujących wektorów kolumnowych:
𝑼 = col(𝑈
1
, 𝑈
2
, … , 𝑈
𝑏
),
(1)
𝑰 = col(𝐼
1
, 𝐼
2
, … , 𝐼
𝑏
),
𝑽 = col(𝑉
1
, 𝑉
2
, … , 𝑉
𝑛
),
gdzie 𝑽, to wektor potencjałów węzłowych czyli potencjałów 𝑛 węzłów w odniesieniu do potencjału węzła odniesienia
(czyli węzła o numerze zero).
Prądy i napięcia we wszystkich gałęziach obwodu złożonego z 𝑏 gałęzi (na przykład o postaci jak na rysunku 1b) oraz
𝑛 węzłów niezależnych i węzła odniesienia, podlegają prawom Kirchhoffa, które można przedstawić za pomocą nastę-
pujących zapisów macierzowych
𝑨 ⋅ 𝑰 = 𝟎
(2)
𝑼 = 𝑨
𝑻
⋅ 𝑽
gdzie 𝑨
𝑻
to transponowana macierz incydencji.
Do pełnego opisu obwodu brakuje równań opisujących zachowanie się elementów; dla oporników są to równania prawa
Ohma, dla niezależnych idealnych źródeł napięciowych równania postaci 𝑈
𝐸
= 𝐸, itd. Takich równań powinno być łącz-
nie 𝑏. W zapisie macierzowym można im nadać formę
𝑴 ⋅ 𝑰 + 𝑵 ⋅ 𝑼 = 𝑬𝒙𝒄,
(3)
gdzie 𝑴 = [𝑚
𝑖𝑗
]
𝑏×𝑏
, 𝑵 = [𝑛
𝑖𝑗
]
𝑏×𝑏
, 𝑬𝒙𝒄 = col(𝑒𝑥
1
, 𝑒𝑥
2
, … , 𝑒𝑥
𝑏
).
Na przykład, gdy element o numerze 𝑘 jest opornikiem o oporze 𝑅
𝑘
, to 𝑚
𝑘𝑘
= 𝑅
𝑘
, 𝑛
𝑘𝑘
= −1, 𝑒𝑥
𝑘
= 0, co odpo-
wiada (pod warunkiem, że 𝑚
𝑘𝑗
= 𝑛
𝑘𝑗
= 0 dla 𝑗 ≠ 𝑘) zapisowi prawa Ohma dla tego elementu w postaci
𝑅
𝑘
⋅ 𝐼
𝑘
− 𝑈
𝑘
= 0.
Podobnie, gdy element o numerze 𝑙 jest źródłem prądowym o SPM 𝐽
𝑙
, to 𝑚
𝑙𝑙
= 1, 𝑛
𝑙𝑙
= 0, 𝑒𝑥
𝑘
= 𝐽
𝑙
, co odpowiada
(pod warunkiem, że 𝑚
𝑙𝑗
= 𝑛
𝑙𝑗
= 0 dla 𝑗 ≠ 𝑙) zapisowi definicyjnego równania dla tego elementu w postaci
1 ⋅ 𝐼
𝑙
− 0 ⋅ 𝑈
𝑙
= 𝐽
𝑙
(założono, że strzałki 𝐼
𝑙
oraz 𝐽
𝑙
są zgodne).
Powyższe równania macierzowe dają się zapisać w postaci
[
𝑨
𝟎
𝟎
𝟎
𝑨
𝑇
−𝟏
𝑴 𝟎
𝑵
] ⋅ [
𝑰
𝑽
𝑼
] = [
𝟎
𝟎
𝑬𝒙𝒄
] .
(4)
Macierz kwadratowa z ostatniego równania nosi nazwę macierzy tableau.
Powyższe równanie może być łatwo przekształcone do postaci
[𝑨
𝟎
𝑴 𝑵⋅𝑨
𝑇
] ⋅ [𝑰
𝑽
] = [𝟎
𝑬𝒙𝒄
] .
(5)
Równanie (5) jest macierzowym zapisem układu (𝑏+𝑛) równań liniowych z (𝑏+𝑛) niewiadomymi (na marginesie
zauważmy, że, przy naszych oznaczeniach, 𝑏+𝑛 byłoby liczbą niezależnych cykli w napięciowym grafie obwodu, co ozna-
cza, że ten układ równań jest równoliczny z układem równań, jaki uzyskalibyśmy stosując metodę prądów obwodowych
(w szczególności oczkowych)).
Rozwiązaniem układu równań (5) są potencjały wszystkich (niezależnych) 𝑛 węzłów obwodu oraz prądy dwójników
obwodu. Przy oznaczeniach 𝑺 = [
𝟎
𝑬𝒙𝒄
], 𝑿 = [𝑰
𝑽
], 𝑇 = [𝑨 𝟎
𝑴 𝑵⋅𝑨
𝑇
] rozwiązanie 𝑿 możemy uzyskać jedną instrukcją MA-
TLAB-a:
X=S\T
(6)
Bilans mocy w obwodzie
W poprzednich ćwiczeniach stwierdziliśmy, że bilans mocy w obwodzie układamy przyjmując, że iloczyn prądu i napięcia przy przeciwnej orien-
tacji napięcia na dwójniku i płynącego przez ten dwójnik prądu oznacza moc traconą na tym dwójniku; zawsze na rezystorach moc tak liczona jest
dodatnia, zaś bardzo często liczona w ten sam sposób moc na źródłach jest ujemna; ujemną moc traconą należy interpretować (po zmianie znaku)
jako dodatnią moc dostarczoną.
Spójrzmy na zagadnienie bilansu najpierw ogólnie, pamiętając, że dla prądów i napięć gałęzi mieliśmy zależności
(patrz wzory (2)):
𝑨 ⋅ 𝑰 = 𝟎, 𝑼 = 𝑨
𝑻
⋅ 𝑽.
C. Stefański, Materiały do ćw 6 Lab. TOiS
4/12
Ostatnia modyfikacja: maj 2016
Wykorzystamy je w następującym rachunku:
∑𝑃 = 𝑼
𝑇
⋅ 𝑰 = (𝑨
𝑻
⋅ 𝑽)
𝑇
⋅ 𝑰 = (𝑽
𝑻
⋅ 𝑨) ⋅ 𝑰 = 𝑽
𝑻
⋅ (𝑨 ⋅ 𝑰) = 𝑽
𝑻
⋅ 𝟎 = 0 .
(7)
Oznacza to, że suma mocy traconych na wszystkich elementach w obwodzie o 𝑏 dwójnikach wynosi zero (jak zresztą
wynika z twierdzenia Tellegena).
Bilans mocy winien być oczywiście zerowy (suma po wszystkich elementach mocy traconych na elementach = zero):
"
∑𝑃
dwójników
tracona
= 0.
(8)
Koncepcja źródeł zastępczych Thevenina i Nortona
Jak pamiętamy z wcześniejszego opracowania, jeżeli nasze potrzeby ograniczają się do znalezienia prądu (napięcia, mocy) tylko w kilku zgrupo-
wanych elementach (w skrajnym przypadku w jednym elemencie) obwodu, to całą pozostałą część struktury możemy zastąpić tzw. źródłem Thevenina
lub źródłem Nortona.
Przyjmując, tym razem, że wybranym ważnym elementem jest idealne źródło prądowe 𝐽
5
możemy,
zgodnie z powyższą koncepcją, wieloelementowy obwód zredukować do jednej z przedstawionych na rysunku 3a,b
po-
staci.
)
b
i
k
J
5
5
R
T
E
T
U
I
)
a
(y)
(x)
Obwód liniowy
zastępowany
źródłem
Thevenina
lub Nortona
q
-t
y
w
ar
ia
n
t
o
b
ci
ąż
en
ia
o
b
w
o
d
u
, q
=1
lu
b
2
i
k
J
5
5
R
N
J
N
U
I
)
c
q
U
q
I
Rys. 3. Obwód zastępczy Thevenina a), obwód zastępczy Nortona b), schemat ogólny do wyznaczania parame-
trów obwodu zastępczego c).
Poza parametrami idealnego źródła prądowego 𝐽
5
, parametry pozostałych elementów z rysunków 3a,b muszą zo-
stać dopiero obliczone na tej podstawie, że na źródłe 𝐽
5
winien wystąpić taki sam spadek napięcia 𝑈, jaki występował
na tym elemencie w ramach obwodu z rysunku 1b.
Wprowadzone przez Thevenina oraz Nortona definicje dają nam następujące zapisy, które odnoszą się do ry-
sunku 1b i rysunków 3a,b:
𝐸
𝑇
= 𝑈
5
|
𝐽
5
=0
, 𝐽
𝑁
= 𝐼
5
|
𝐽
5
zastąpione zwarciem
, 𝑅
𝑇
= 𝑅
𝑁
=
𝐸
𝑇
𝐽
𝑁
.
(9a)
Inaczej sprawę ujmując, napięcie Thevenina otrzymujemy wtedy, kiedy w obwodzie z rysunku 1b wyzerujemy ele-
ment 𝐽
5
, co odpowiada rozwarciu zacisków obwodu zastępowanego źródłem Thevenina. W sytuacji analizy metodą ta-
bleau wymaga to tylko zmiany wartości elementu 𝐽
5
bez zmiany kształtu równania (5).
Prąd Nortona otrzymamy na podstawie rozwiązania równań (5) po zastąpieniu prądowego źródła 𝐽
5
źródłem napię-
ciowym 𝐸
5
= 0. Co prawda nie wymaga to tym razem od nas przekonstruowania równań typu (5) opisujących obwód
(bowiem algorytm uczyni to automatycznie), ale dane wejściowe trzeba będzie skorygować.
Na szczęście, jak już wiemy, do wyznaczenia parametrów źródeł Thevenina i Nortona można podejść bardziej jed-
nolicie, mianowicie parametry te można obliczyć z ogólniejszych zależności, przy oznaczeniach z rysunku 3c:
𝐸
𝑇
=
𝐼
2
⋅ 𝑈
1
− 𝐼
1
⋅ 𝑈
2
𝐼
2
− 𝐼
1
, 𝐽
𝑁
=
𝑈
2
⋅ 𝐼
1
− 𝑈
1
⋅ 𝐼
2
𝑈
2
− 𝑈
1
, 𝑅
𝑇
= 𝑅
𝑁
=
𝑈
1
− 𝑈
2
𝐼
2
− 𝐼
1
(9b)
Wynikiem analizy tableau naszego obwodu są prądy i napięcia elementowe (a także potencjały węzłowe), zatem
wzór (9b) może być wykorzystany bezpośrednio. Na przykład można przyjąć, że ( 𝑈
1
, 𝐼
1
) są odpowiednio napięciem
i prądem źródła 𝐽
5
w obwodzie bez modyfikacji, zaś ( 𝑈
2
, 𝐼
2
) to odpowiednio napięcie i prąd źródła 𝐽
5
, gdy wartość 𝐽
5
=
0 :
𝐸
𝑇
= 𝑈
2
, 𝐽
𝑁
=
𝑈
2
⋅ 𝐼
1
𝑈
2
− 𝑈
1
, 𝑅
𝑇
= 𝑅
𝑁
=
𝑈
2
− 𝑈
1
𝐼
1
.
(9c)
Znowu do wyznaczenia źródeł Thevenina lub Nortona, gdy korzystamy z wzorów (9b, c) wystarczy odpowiednio
manipulować „obciążeniem”, a nie strukturą całego obwodu.
Policzmy jeszcze moce w obwodach 3a-b, starając się trzymać w konwencji, że korzystamy z wyliczonych prądów
i napięć elementowych na dopełniającym źródło zastępcze dwójniku 𝐽
5
.
Moce (liczone jako tracone (!!!)) w obwodzie z rysunku 3a
C. Stefański, Materiały do ćw 6 Lab. TOiS
5/12
Ostatnia modyfikacja: maj 2016
𝑃
rezystorów
tracona
= 𝐽
5
2
⋅ 𝑅
𝑇
,
(13)
𝑃
źródeł
tracona
= −𝐽
5
⋅ 𝐸
𝑇
+ 𝐽
5
⋅ 𝑈
5
.
Moce (liczone jako tracone (!!!)) w obwodzie z rysunku 3b
𝑃
rezystorów
tracona
=
𝑈
5
2
𝑅
𝑁
,
(14)
𝑃
źródeł
tracona
= −𝐽
𝑁
⋅ 𝑈
5
+ 𝐽
5
⋅ 𝑈
5
.
Warto zastanowić się, czy (i dlaczego) suma mocy traconych na wszystkich, oprócz idealnego źródła prądowego 𝐽
5
,
elementach obwodu 1b jest równa sumie mocy traconych na wszystkich, oprócz idealnego źródła prądowego 𝐽
5
, ele-
mentach obwodu z rysunku 3a, jak i 3b.
3. Wspomagane programem MATLAB rozwiązywanie układu równań liniowych
Przypomnimy teraz krótko to, co zostało bardziej szczegółowo omówione w instrukcji do poprzedniego ćwiczenia, a co dotyczyło rozwiązywania
układu równań liniowych wspomaganego narzędziami typu MATLAB (lub mu podobnymi, jak OCTAVE, czy SciLab). Układ równań można zapisywać w
formie „klamrowej” (15):
{
𝑎
11
𝑥
1
+ 𝑎
12
𝑥
2
+ ⋯ 𝑎
1𝑛
𝑥
𝑛
= 𝑏
1
𝑎
21
𝑥
1
+ 𝑎
22
𝑥
2
+ ⋯ 𝑎
2𝑛
𝑥
𝑛
= 𝑏
2
… … … … … … … … … … .
𝑎
𝑛1
𝑥
1
+ 𝑎
𝑛2
𝑥
2
+ ⋯ 𝑎
𝑛𝑛
𝑥
𝑛
= 𝑏
𝑛
(15)
lub w postaci macierzowej
[
𝑎
11
𝑎
12
𝑎
13
… . 𝑎
1𝑁
𝑎
21
𝑎
22
𝑎
23
… . 𝑎
2𝑁
𝑎
31
𝑎
32
𝑎
33
… 𝑎
3𝑁
…
…
…
…
…
𝑎
𝑁1
𝑎
𝑁2
𝑎
𝑁3
… 𝑎
𝑁𝑁
]
⋅
[
𝑥
1
𝑥
2
𝑥
3
…
𝑥
𝑁
]
=
[
𝑏
1
𝑏
2
𝑏
3
…
𝑏
𝑁
]
.
(16)
Skrócony zapis (przy notacji bez nawiasów kwadratowych) równania (16) wygląda następująco:
𝑨 ∙ 𝒙 = 𝒃
(17)
W tym przypadku rozwiązanie jest w naszym zasięgu nawet bez wspomagania komputerowego, przy czym często korzystamy z przejrzystych
zapisów metody Cramera
𝒙
𝟏
=
𝜟
𝟏
𝜟
, 𝒙
𝟐
=
𝜟
𝟐
𝜟
, 𝒙
𝟑
=
𝜟
𝟑
𝜟
, 𝒙
𝟒
=
𝜟
𝟒
𝜟
, 𝜟 ≠ 𝟎 .
(18)
gdzie Δ jest wyznacznikiem macierzy 𝑨, zaś każdy z Δ
𝑛
to wyznacznik macierzy powstałej z 𝑨 przez zastąpienie w niej wektorem 𝒃 kolumny o numerze
„𝑛-tym”.
Współcześnie, spośród wielu opracowanych metod numerycznego rozwiązywania układów równań liniowych, duże znaczenie mają różne wa-
rianty metody eliminacji Gaussa. Algorytm tej metody dla równania (16) jest następujący (eliminacja niewiadomych w wyrażeniu (16) przebiega
następująco):
W pierwszym kroku od wierszy 2, 3, … . , 𝑁 odejmujemy wiersz pierwszy pomnożony przez taki współczynnik, który zapewnia jako
wynik odejmowania nowy wiersz nie zawierający już niewiadomej 𝑥
1
.
Po pierwszym kroku (krok 𝑠 = 1; ang.: step) otrzymujemy
[
𝑎
11
𝑎
12
𝑎
13
… 𝑎
1𝑁
0 𝑎
22
(1)
𝑎
23
(1)
… 𝑎
2𝑁
(1)
0 𝑎
32
(1)
𝑎
33
(1)
… 𝑎
3𝑁
(1)
… …
…
… …
0 𝑎
𝑁2
(1)
𝑎
𝑁3
(1)
… 𝑎
𝑁𝑁
(1)
]
⋅
[
𝑥
1
𝑥
2
𝑥
3
…
𝑥
𝑁
]
=
[
𝑏
1
𝑏
2
(1)
𝑏
3
(1)
…
𝑏
𝑁
(1)
]
(19a)
gdzie
𝑎
𝑖𝑘
(1)
= 𝑎
𝑖𝑘
− 𝑙
𝑖1
⋅ 𝑎
1𝑘
,
𝑏
𝑖
(1)
= 𝑏
𝑖
− 𝑙
𝑖1
⋅ 𝑏
1
, przy 𝑙
𝑖1
=
𝑎
𝑖1
𝑎
11
oraz 𝑖, 𝑘 = 2, 3, … , 𝑁
W drugim kroku od wierszy 3, 4, … . , 𝑁 odejmujemy wiersz drugi pomnożony przez taki współczynnik, który zapewnia jako wynik
odejmowania nowy wiersz nie zawierający już niewiadomej 𝑥
2
.
Po drugim kroku (krok 𝑠 = 2) otrzymujemy
[
𝑎
11
𝑎
12
𝑎
13
… 𝑎
1𝑁
0 𝑎
22
(1)
𝑎
23
(1)
… 𝑎
2𝑁
(1)
0 0
𝑎
33
(2)
… 𝑎
3𝑁
(2)
… …
…
… …
0 0
𝑎
𝑁3
(2)
… 𝑎
𝑁𝑁
(2)
]
⋅
[
𝑥
1
𝑥
2
𝑥
3
…
𝑥
𝑁
]
=
[
𝑏
1
𝑏
2
(1)
𝑏
3
(2)
…
𝑏
𝑁
(2)
]
(19b)
gdzie
𝑎
𝑖𝑘
(2)
= 𝑎
𝑖𝑘
(1)
− 𝑙
𝑖2
⋅ 𝑎
2𝑘
(1)
,
𝑏
𝑖
(2)
= 𝑏
𝑖
(1)
− 𝑙
𝑖2
⋅ 𝑏
2
(1)
, przy 𝑙
𝑖2
=
𝑎
𝑖2
(1)
𝑎
22
(1)
oraz 𝑖, 𝑘 = 3, 4, … , 𝑁
C. Stefański, Materiały do ćw 6 Lab. TOiS
6/12
Ostatnia modyfikacja: maj 2016
Powyższe postępowanie stosujemy do otrzymanego w drugim kroku równania (20b) i kolejno w następnych krokach do następnych
równań ze zwiększającą się ilością zer pod przekątną główną przetwarzanej macierzy; po (N-1) krokach otrzymujemy
[
𝑎
11
𝑎
12
𝑎
13
… 𝑎
1𝑁
0
𝑎
22
(1)
𝑎
23
(1)
… 𝑎
2𝑁
(1)
0
0
𝑎
33
(2)
… 𝑎
3𝑁
(2)
…
…
…
… …
0
0
0
… 𝑎
𝑁𝑁
(𝑁−1)
]
⋅
[
𝑥
1
𝑥
2
𝑥
3
…
𝑥
𝑁
]
=
[
𝑏
1
𝑏
2
(1)
𝑏
3
(2)
…
𝑏
𝑛
(𝑁−1)
]
(20)
gdzie
𝑎
𝑖𝑘
(𝑠)
= 𝑎
𝑖𝑘
(𝑠-1)
− 𝑙
𝑖𝑠
⋅ 𝑎
𝑠𝑘
(𝑠-1)
, 𝑏
𝑖
(𝑠)
= 𝑏
𝑖
(𝑠-1)
− 𝑙
𝑖𝑠
⋅ 𝑏
𝑠
(𝑠-1)
,
przy 𝑙
𝑖𝑠
=
𝑎
𝑖𝑠
(𝑠-1)
𝑎
𝑠𝑠
(𝑠-1)
, 𝑠 = 1, 2, … , 𝑁-1, 𝑖 = 𝑠+1, 𝑠+2, … , 𝑁, 𝑘 = 𝑖-1, 𝑖, … 𝑁 (przy czym 𝑎
𝑝𝑞
(0)
= 𝑎
𝑝𝑞
).
Przykład obliczeń według tego algorytmu był podany w instrukcji do poprzedniego ćwiczenia – warto z nim ponownie zapoznać się.
Poniżej przypominamy też pewien wariant metody Gaussa noszący nazwę metody LU. Polega ona na przyjęciu lewostronnej macierzy ze wzoru
(20) jako macierzy 𝑼 (trójkątnej górnej) i utworzeniu dodatkowej macierzy 𝑳 (trójkątnej dolnej) zbudowanej ze współczynników 𝑙
𝑖𝑠
. Przykładowo dla
przypadku macierzy o wymiarach 4x4 (𝑁 = 4) otrzymujemy
𝑼 =
[
𝑎
11
𝑎
12
𝑎
13
𝑎
14
0
𝑎
22
(1)
𝑎
23
(1)
𝑎
24
(1)
0
0
𝑎
33
(2)
𝑎
34
(2)
0
0
0
𝑎
44
(3)
]
, 𝑳 =
[
1
0
0
0
𝑙
21
1
0
0
𝑙
31
𝑙
32
1
0
𝑙
41
𝑙
42
𝑙
43
1]
(21)
Macierze te, dla 𝑁 ≥ 2, mają następujące pożyteczne właściwości
𝑳 ∙ 𝑼 = 𝑨, |𝑳| = 1, |𝑨| = |𝑼|
(22)
Wzory (23) i (24) wskazują na korzystniejszy od tradycyjnego sposób obliczania wyznacznika macierzy 𝑨; dodajmy, sposób wykorzystywany w
operacji det(A) MATLAB-a.
Poszukiwanie wektora rozwiązań 𝒙 możemy teraz przedstawić w postaci
𝑳 ∙ 𝑼 ∙ 𝒙 = 𝒃 ,
(23)
𝒚 = 𝑳
−𝟏
∙ 𝒃 ,
𝒙 = 𝑼
−𝟏
∙ 𝒚 ,
gdzie 𝒚 to pomocniczy wektor kolumnowy w procedurze (23).
Eliminację Gaussa prowadzącą do otrzymania macierzy LU można także zastosować nie do oryginalnego układu równań (15), tylko do układu
równań z przestawionymi wierszami. Najbardziej oczywista potrzeba przestawienia wierszy może w ramach algorytmu eliminacji wynikać z potrzeby
ominięcia sytuacji dzielenia przez zero, ale możemy także w ten sposób zwiększać dokładność obliczeń.
Wspomnianą optymalizację dokładności obliczeń w odniesieniu do układu równań (17) otrzymujemy w MATLAB-ie po zastosowaniu do dekom-
pozycji funkcji lu MATLAB-a z argumentami wyjściowymi 𝑳, 𝑼, 𝑷:
[L,U,P] = lu(A).
(24)
Funkcja ta w miejsce macierzy 𝑨 zwraca nam trzy macierze, które możemy podstawić do równania (17). Po podstawieniu, zapis tego równania przyj-
muje postać
𝑷 ∙ 𝑳 ∙ 𝑼 ∙ 𝒙 = 𝒃
.
(25)
Zawierająca informację o przestawieniach wierszy macierz 𝑷 ma następujące właściwości
𝑷 = 𝑷
−𝟏
,
|𝑷| = 1 .
(26)
Stąd rozwiązanie równania (17) możemy w tym przypadku uzyskać na podstawie wyrażenia
𝑳 ∙ 𝑼 ∙ 𝒙 = 𝒃
𝒑
(27)
gdzie
𝒃
𝒑
= 𝑷 ∙ 𝒃.
Sposób zamiany kolejności wierszy daje się odczytać z macierzy 𝑷. Ściśle rzecz biorąc, niezerowy wyraz 𝑝
𝑖𝑗
= 1 tej macierzy informuje nas, że
nowym 𝑖-tym równaniem w zapisie (15) ma być stare 𝑗-te równanie tego zapisu.
Realizując poleceniem lu(A)opisaną powyżej procedurę stosujemy optymalizowaną dekompozycję [L, U, P] w sposób jawny, ale oprócz
tego dekompozycja taka jest automatycznie wykonywana niejawnie w ramach następujących poleceń Matlaba:
w przypadku równania postaci (16), gdy je rozwiązujemy poleceniem
x=A\b
(28)
w przypadku obliczania wyznacznika macierzy A
DA = det(A)
(29)
5. Pliki MATLABA do wspomagania analizy obwodów
Przyjęto podobny do stosowanego w programach typu SPICE sposób zadawania struktury i wartości elementów
obwodu czyli strukturę danych o obwodzie. Sposób ten w literaturze bywa nazywany tablicą połączeń.
C. Stefański, Materiały do ćw 6 Lab. TOiS
7/12
Ostatnia modyfikacja: maj 2016
Zakładamy, że wiersz danych ma postać tekstową o następującej strukturze:
ReNe nv1 nv2 val
gdzie: Re to rodzaj elementu (na razie jedno z {R, G, E, J}), Ne to nr elementu (wymaga się ciągłej numeracji elementów
poczynając od 1); nv1 i nv2 to pierwszy i drugi węzeł, między które załączono element (wymaga się ciągłej numeracji
węzłów poczynając od zera (węzeł odniesienia)), wreszcie val oznacza wartość liczbową elementu; dopuszcza się krot-
ności {p,n,u,m,k,M,G,T}.
Na jeden dwójnik przypada jeden wiersz tablicy połączeń. Przy pisaniu algorytmów założono, że nv1 oznacza węzeł
o umownie wyższym potencjale (umowny kierunek prądu od nv1 ku nv2).
Aby nie komplikować oprogramowania przyjęto, że liczba węzłów i gałęzi jest podawana na początku danych (pro-
fesjonalne programy takiego wymagania nie mają).
% cramer.m
function [x]=cramer(A,b)
%Metoda Cramera rozwiązywania układów równań liniowych
%Zapis układu równań A*x=y
N=rank(A);
for k=1:N
K=A(:,k); % Kolumna o numerze k w macierzy A na starcie
A(:,k)=b; % Macierze A
1
, A
2
,…,A
m
DK(k)=det(A); % Wyznaczniki D
1
, D
2
,…,D
m
A(:,k)=K; %Powrót do startowej postaci macierzy A
end
D=det(A); % Wyznacznik główny
x=DK./D; % Tablicowe dzielenie wyznaczników
%Koniec cramer.m
% dekompozycja.m
function [LA,UA]=dekompozycja(A)
% Dekompozycja (rozkład) macierzy A na dwie macierze
%„trójkątną-dolną” LA i „trójkątną- górną” UA
N=rank(A); % Określa liczbę niezależnych równań
el=eye(N); % Generuje macierz z "jedynkami" na głównej przekątnej
a=zeros(N); % Generuje macierz wypełnioną zerami
a=A; % Przepisuje elementy macierzy
for s=1:(N-1) % Litera s oznacza numer kroku w metodzie Gaussa
% Litery i oraz k oznaczają numery wiersza i kolumny
n=s+1;
for i=n:N
el(i,s)=a(i,s)/a(s,s)
for k=s:N
a(i,k)=a(i,k)-el(i,s)*a(s,k);
end
end
end
LA=el; % Macierz trójkątna dolna
UA=a; % Macierz trójkątna górna
% Koniec dekompozycja.m
% DaneOobwodzie.m
%Do Cw4
%***Program wspomagający wpisywanie tabeli połączeń elementów***
%
% Nową tabelę połączeń realizujemy poprzez edycję niniejszej
% tabeli połączeń
%
% Dane o obwodzie
TablicaPolaczen=[
"* to jest proba-test";
"E1 1 0 1";
"R2 2 1 1";
"* to jest proba-test";
"R3 2 0 1";
"*koniec";
"R4 3 2 1";
"R5 3 0 1";
"*koniec";
C. Stefański, Materiały do ćw 6 Lab. TOiS
8/12
Ostatnia modyfikacja: maj 2016
"J6 3 4 1";
"R7 4 0 1";
"R8 4 5 1";
"R9 5 0 1";
"R10 5 6 1";
"E11 6 0 1";
"koniec"];
%
% UzupelnijMacierzIncydencji.m
function UzupelnijMacierzIncydencji(A,d)
global A;
if (d(2)>0)
A(d(2),d(1))=1;
endif
if (d(3)>0)
A(d(3),d(1))=-1;
endif
endfunction
% UzupelnijWpisyEWmacierzachNiExc.m
function UzupelnijWpisyEWmacierzachNiExc(N,Exc,d)
global N;
global Exc;
if (d(1)>0)
N(d(1),d(1))=1;
Exc(d(1),1)=d(4);
else
print("4E. Blad danych\n");
close;
endif
endfunction
% UzupelnijWpisyGWmacierzachMiN.m
function UzupelnijWpisyGWmacierzachMiN(M,N,d)
global M;
global N;
if (d(1)>0)
M(d(1),d(1))=-1;
N(d(1),d(1))=d(4);
else
print("4G. Blad danych\n");
close;
endif
endfunction
% UzupelnijWpisyJWmacierzachMiExc.m
function UzupelnijWpisyJWmacierzachMiExc(M,Exc,d)
global M;
global Exc;
if (d(1)>0)
M(d(1),d(1))=1;
Exc(d(1),1)=d(4);
else
print("4J. Blad danych\n");
close;
endif
endfunction
% UzupelnijWpisyRWmacierzachMiN.m
function UzupelnijWpisyRWmacierzachMiN(M,N,d)
global M;
C. Stefański, Materiały do ćw 6 Lab. TOiS
9/12
Ostatnia modyfikacja: maj 2016
global N;
if (d(1)>0)
M(d(1),d(1))=d(4);
N(d(1),d(1))=-1;
else
print("4R. Blad danych\n");
close;
endif
endfunction
% EJRG.m
% Pojedynczy wiersz z danymi oznaczamy WDanych
%
% Zakładamy, że wiersz danych ma postać tekstową o następującej strukturze:
%
% ReNe nv1 nv2 val
% gdzie
% Re to rodzaj elementu (na razie jedno z {R, G, E, J})
% Ne to nr elementu; wymaga się ciągłej numeracji elementów poczynając od 1;
% nv1 i nv2 to pierwszy i drugi węzeł, między które załączono element;
% wymaga się ciągłej numeracji węzłów poczynając od zera (węzeł odniesienia)
% val oznacza wartość liczbową elementu; dopuszcza się krotności {p,n,u,m,k,M,G,T}
global A;
global M;
global N;
global Exc;
LE=input("Liczba elementow wynosi (numeracja od 1 do LE): LE=");
NN=input("Wezel o najwyzszym numerze to (numeracja od 0 do NN): NN=");
A=zeros(NN,LE);
M=zeros(LE,LE);
N=M;
Exc=zeros(LE,1);
T=zeros(LE+NN,LE+NN);
S=zeros(LE+NN,1);
RElem=0;
i=1;
k=0;
DaneOobwodzie;
while (abs(RElem-5)>0.1)
if (i>LE)
break;
else
k=k+1;
endif;
% printf("%d %d \n",i,k);
WDanych=TablicaPolaczen(k,:);
zp=WDanych;
X = strsplit(zp);
if (zp(1,1)=="R")
RElem=1;
elseif (zp(1,1)=="G")
RElem=2;
elseif (zp(1,1)=="E")
RElem=3;
elseif (zp(1,1)=="J")
RElem=4;
elseif (zp(1,1:6)=="koniec")
RElem=5;
elseif (zp(1,1)=="*")
RElem=6;
else
printf ("1. Nierozpoznany element\n");
close;
endif;
C. Stefański, Materiały do ćw 6 Lab. TOiS
10/12
Ostatnia modyfikacja: maj 2016
if (and(RElem==5,i<LE))
printf("2a. Blad danych; nacisnij Ctrl+C\n");
keyboard;
endif;
if(and(size(X)<4, RElem<5))
printf ("2b. Blad danych\n");
close;
endif;
if (RElem<6)
b=strjoin(X(1,1));
a=b;
a(1)=[];
d(1)=str2num(a);
d(2)=str2num(strjoin(X(1,2)));
d(3)=str2num(strjoin(X(1,3)));
c = strrep(strrep (X(1,4), "m", "e-3"),"k","e3");
c = strrep(strrep (c, "u", "e-6"),"M","e6");
c = strrep(strrep (c, "n", "e-9"),"G","e9");
c = strrep(strrep (c, "p", "e-12"),"G","e12");
d(4)=str2num(strjoin(c));
i=i+1;
if(or(eq(d(2),d(3)), gt(d(2),NN), gt(d(3),NN), gt(d(1),LE)))
printf ("3. Blad danych\n");
close;
endif;
if (RElem==1)
UzupelnijWpisyRWmacierzachMiN(M,N,d);
endif;
if (RElem==2)
UzupelnijWpisyGWmacierzachMiN(M,N,d);
endif;
if (RElem==3)
UzupelnijWpisyEWmacierzachNiExc(N,Exc,d);
endif;
if (RElem==4)
UzupelnijWpisyJWmacierzachMiExc(M,Exc,d);
endif;
if (RElem <= 4)
UzupelnijMacierzIncydencji(A,d(1:3));
endif;
endif;
endwhile;
T(1:NN,1:LE)=A;
T(NN+1:NN+LE,1:LE)=M;
T(NN+1:NN+LE,LE+1:LE+NN)=N*A';
S(NN+1:NN+LE,1)=Exc;
x=T\S;
I=x(1:LE,1);
Vn=x(LE+1:LE+NN,1);
U=A'*Vn;
I',U',Vn'
%
6. Symulacja obwodu elektrycznego za pomocą programu SPICE (Micro̵CAP)
Z wielu wariantów SPICE’a w laboratorium jest do dyspozycji program Micro̵CAP w wersji 10 lub 11 (dalej będziemy
na oznaczenie tego narzędzia używać skrótu MC).
C. Stefański, Materiały do ćw 6 Lab. TOiS
11/12
Ostatnia modyfikacja: maj 2016
W edytorze programu MC tworzymy wirtualny obwód złożony: ze źródeł niezależnych, źródeł sterowanych, rezy-
storów, cewek indukcyjnych i kondensatorów i w miarę potrzeby z innych elementów pochodzących z banku elementów
SPICE’a. SPICE – w zewnętrznym oglądzie – raczej symuluje właściwości obwodu niż wspomaga poszukiwanie rozwiązań
równań opisujących taki obwód, choć – w ukryciu przed użytkownikiem – naprawdę te równania automatycznie tworzy
i jedno- lub wielokrotnie te równania rozwią-
zuje.
Obwód prądu stałego z rys. 1b, wprowa-
dzony do MC, wygląda jak na rys. 4.
Edytor MC nie pozwala w wygodny sposób
określać kierunki odniesienia dla przepływu prą-
dów w gałęziach oraz nadawać numery porząd-
kowe węzłów. We wprowadzonym do edytora
programu obwodzie można poprzez kliknięcie
myszką w obrębie wybranego elementu zmienić
w otwierającym się wtedy okienku wartość licz-
bową parametru (np. E1, J5, R2, R3). Do równoczesnego i sprawnego nadawania wartości wielu elementom zdefinio-
wano stałą symboliczną N (linia .define N 30). Od wartości tej stałej uzależnione są tu SEM źródła E1, SPM źródła
J5, oporności oporników R2, R3 i R4.
Jeżeli z kilku wariantów badania obwodu DC (prądu stałego) wybierzemy analizę Dynamic DC (AnalysisDy-
namic DC…) otrzymamy napięcia węzłowe (V1, V2), czyli potencjały każdego z ponumerowanych węzłów względem
napięcia odniesienia (masa, ziemia) oraz prądy gałęziowe (I1, I2, …, I5). Na ekran wartości tych napięć i prądów można
wywołać wybierając odpowiednie ikonki z paska narzędziowego programu
1
.
Oprócz dotychczas pokazanych sposobów wyznaczania parametrów źródeł zastępczych Thevenina lub Nortona
można też posłużyć się następującymi podejściami, sugerowanymi rysunkiem 5, dość wygodnymi, gdy korzystamy z pro-
gramu MC. Mianowicie analiza Dynamic DC przy 𝐽
𝑡𝑒𝑠𝑡
=0 pozwala wyznaczyć 𝐸
𝑇
bezpośrednio, zaś przy 𝐽
𝑡𝑒𝑠𝑡
=-1 A licz-
bowa wartość napięcia 𝑈
𝐽
𝑡𝑒𝑠𝑡
różni się liczbowo o 𝑅
𝑇
od liczbowej wartości 𝐸
𝑇
; podobnie analiza Dynamic DC przy
𝐸
𝑡𝑒𝑠𝑡
= 0 pozwala wyznaczyć 𝐽
𝑁
bezpośrednio, zaś przy 𝐸
𝑡𝑒𝑠𝑡
= −1 V liczbowa wartość prądu 𝐼
𝐸
𝑡𝑒𝑠𝑡
różni się liczbowo
o 𝐺
𝑁
od liczbowej wartości 𝐽
𝑁
.
Rys. 5. Alternatywne sposoby wyznaczania parametrów źródeł zastępczych
Kontynuując myśl sugerowaną rysunkiem 5 możemy pobudzić badany obwód (ważne, by był to układ zbudowany
tylko z elementów E, J, R, ZS, WO) ze żródła 𝐽
𝑡𝑒𝑠𝑡
o przebiegu okresowym w postaci fali prostokątnej unipolarnej, o mi-
nimalnej wartości -1 A (minus jeden amper) i wartości maksymalnej 0 A. W przedziałach czasu, w których przebieg 𝐽
𝑡𝑒𝑠𝑡
przyjmuje wartość zero, napięcie na źródle 𝐽
𝑡𝑒𝑠𝑡
wyniesie 𝐸
𝑇
(będzie szukanym napięciem rozwarcia), zaś zmiana war-
tości napięcia na źródle 𝐽
𝑡𝑒𝑠𝑡
przy zmianie prądu na wartość minimalną (-1 A) liczbowo będzie równa rezystancji Theve-
nina. Badanie to w Micro∙Cap-ie przeprowadzić można stosując analizę TRANSIENT. Odpowiednie ilustracje rysunku 6
przybliżają ten pomysł dla przypadku, gdy szukamy źródła Thevenina dla układu bez R4 z punktu widzenia zacisków (2)
i (0).
1
W MC11 ikonki te wyglądają, jak pokazano tutaj
.
(a)
U
Jtest
=E
T
-R
T
·J
test
R
T
·J
test
J
test
R
T
E
T
U
Jtest
=-R
T
·
J
test
J
test
=0
U
Jtest
=E
T
J
test
=-1A
U
Jtest
=R
T
·1A
I
Etest
=J
N
-G
N
·E
test
G
N
·E
test
E
test
J
N
I
Etest
=-G
N
·
E
test
E
test
=0
I
Etest
=J
N
E
test
=-1V
I
Etest
=G
N
·1V
(a)
(b)
(b)
G
N
E1
1.11*N
R2
2.22*N
R3
3.33/N
J5
DC 5.55/N
Jtest
.define N 30
1
2
Rys. 4. Schemat obwodu z rys. 1b generowany przez program
Micro̵CAP
E1
1.11*N
R2
2.22*N
R3
3.33/N
R4
4.44/N
J5
DC 5.55/N
.define N 30
1
2
C. Stefański, Materiały do ćw 6 Lab. TOiS
12/12
Ostatnia modyfikacja: maj 2016
Rys. 6. Obwód, okienko nastaw analizy TRANSIENT i jej rezultaty dla wyznaczenia 𝑹
𝑻
i 𝑬
𝑻
Na kolejnym rysunku porównano – uzyskane w programie Micro∙Cap – napięcie na R4=4.4/N omów w obwodzie
oryginalnym (po lewej) i w obwodzie z zastępczym źródłem Thevenina (po prawej).
Rys. 7. Napięcie na R4 w badanym obwodzie i w obwodzie ze źródłem Thevenina
Sugeruje się by Czytelnik zaproponował analogiczną procedurę dla znalezienia 𝐺
𝑁
i 𝐽
𝑁
.
Pytania kontrolne
1. Na czym polega tzw. metoda LU rozwiązywania układu równań liniowych?
2. Co to jest macierz incydencji obwodu? Narysuj graf o zredukowanej macierzy incydencji
𝐴 = [
1
0 −1 1
0
−1
1
0 0
0
0 −1
1 0 −1
]
3. Podać szczegóły tworzenia tzw. źródła Thevenina dla wieloelementowego obwodu liniowego.
4. Podać szczegóły tworzenia tzw. źródła Nortona dla wieloelementowego obwodu liniowego.
5. Podać szczegóły metody rozwiązywania układu równań liniowych związaną z nazwiskiem Cramera.
6. W jaki sposób wprowadzamy dane badanego obwodu do programu Micro̵CAP?
7. Na czym polega tzw. dekompozycja LU macierzy współczynników obwodu?
8. Opisz, jaką ma składnię, jak jest interpretowane i do czego można wykorzystać polecenie .define
9. Objaśnij poszczególne równania macierzowe: 𝑨 ⋅ 𝑰 = 𝟎, 𝑼 = 𝑨
𝑻
⋅ 𝑽, 𝑴 ⋅ 𝑰 + 𝑵 ⋅ 𝑼 = 𝑬𝒙𝒄. W szczególności
podaj również rozmiary macierzy i wektorów, gdy wiesz, że dotyczą one obwodu złożonego tylko z dwójników
(w liczbie 𝑏) połączonych w węzłach (węzłów jest 𝑛).
10. Co to jest tablica połączeń? Objaśnij szczegóły wiersza takiej tablicy na przykładzie elementu wybranego spo-
śród {E, J, R, G}.
0.00u
5.00u
10.00u
15.00u
20.00u
25.00u
-1.00
-0.75
-0.50
-0.25
0.00
0.25
I(Jtest) (A)
T (Secs)
0.00u
5.00u
10.00u
15.00u
20.00u
25.00u
30.00m
60.00m
90.00m
120.00m
150.00m
180.00m
v(2) (V)
T (Secs)
Micro-Cap 10 Evaluation Version
ZnalezienieThevenin.cir
110.815308m
320n,34.906822m
E1
1.11*N
R2
2.22*N
R3
3.33/N
R4
4.44/N
J5
DC 5.55/N
Jtest
1.11*N
2.22*N
3.33/N
DC 5.55/N
4.44/N
ET
34.906822m
RT
110.815308m
.define N 30
33.3
34.906822m
19.96099m
33.3
19.96099m
34.906822m