Ćwiczenie 4 – Programowanie Liniowe cd
Opracował: K. Dziedzicki
Rozwiązać problemy programowania liniowego
• Rozwiązując wymaganą liczbę podukładów równań, eliminując rozwiązania niedopuszczalne oraz
wyznaczając rozwiązanie o najlepszej wartości funkcji celu.
• Budując tablicę simplexową a następnie krok po kroku, przechodząc poprzez kolejne rozwiązania
dopuszczalne wyznaczyć rozwiązanie optymalne
Porównać efektywność obydwu podejść.
a – ilość liter w nazwisku
Zadanie 1
dana jest funkcja celu
f(X) = x
1
+ ax
2
i ograniczenia
x
1
≥
1
x
2
≥
2
7/9x
1
+ x
2
≤
10
1/2x
1
- x
2
≤
1
wyznaczyć wartości zmiennych decyzyjne (x
1
, x
2
) dla których funkcja celu osiąga minimum
Zadanie 2
Przedsiębiorstwo produkcyjne wytwarza 2 rodzaje produktów A i B, w ciągu godziny może
wyprodukować 14t wyrobu A lub 7t wyrobu B. Posiadany transport pozwala na przewiezienie w ciągu
godziny do 7t produktu A lub do 12t produktu B. Urządzenia załadowcze pozwalają załadować nie więcej
niż 8t dowolnego z produktów. Przedsiębiorstwo otrzymuje 5000 PLN/t produktu A i a*1000 PLN/t B.
Jaką ilość produktów przedsiębiorstwo powinno produkować aby uzyskać maksymalny dochód. Założyć
że nie istnieje problem sprzedaży produkcji.
Przykład
Elektrownia
składa się z dwóch wydziałów, wydziału produkującego energię elektryczną
i wydziału produkującego parę.
Wyprodukowanie 1t pary wymaga:
2.1
j.p.
środków finansowych,
0.4
węgla brunatnego,
0.045MWh
energii
elektrycznej,
zaś 1MWh
6
j.p.
środków finansowych,
5t
pary
wodnej.
Jedna tona węgla kosztuje 4 j.p., zaś cena sprzedaży 1t pary wodnej wynosi 6j.p, a 1MWh 40j.p.
Określić plan produkcji zapewniający maksymalny zysk, jeśli zasoby środków finansowych wynoszą
600 j.p. a węgla brunatnego 800t.
Oznaczamy
przez
x
1
ilość produkowanej pary a przez x
2
produkcje energii elektrycznej.
Elektrociepłownia powinna wytwarzać parę i energię elektryczną na sprzedaż więc:
x
1
– 5x
2
≥ 0
x
2
– 0.045x
1
≥ 0
zużycie węgla
0.4x
1
≤ 800
wartość zużytych środków finansowych
2.1x
1
+ 6x
2
≤ 600
produkcja nie może być ujemna
x
1
, x
2
≥ 0
Uzyskany efekt finansowy jest różnicą między środkami uzyskanymi ze sprzedaży a kosztami
wytwarzania:
f = 6(x
1
– 5x
2
) + 40(x
2
– 0.045x
1
) – (2.1x
1
+ 4(0.4x
1
) + 6x
2
) =
= 0.5x
1
+ 4x
2
Podsumowując
max { 0.5x
1
+ 4x
2
}
x
1
, x
2
∈
Ω
x
1
– 5x
2
≥ 0
Ω
:
x
2
– 0.045x
1
≥ 0
0.4x
1
≤ 800
2.1x
1
+ 6x
2
≤ 600
x
1
≥ 0
x
2
≥ 0
Rys. 1 Przestrzeń optymalizacji z naniesionymi ograniczeniami
Rozwiązanie 1
lp równania
przecięcie w
spełnia ograniczenia
1
-x
1
+ 5x
2
≤ 0
0.045x
1
-1x
2
≤ 0
x
1
=0.000
x
2
=0.000
tak
2
-x
1
+ 5x
2
≤ 0
0.4x
1
≤ 800
x
1
=2000.000
x
2
=400.000
nie
3
-x
1
+ 5x
2
≤ 0
2.1x
1
+ 6x
2
≤ 600
x
1
=181.818
x
2
=36.364
tak
4
-x
1
+ 5x
2
≤ 0
-x
1
≤ 0
x
1
=0.000
x
2
=0.000
tak
5
-x
1
+ 5x
2
≤ 0
-x
2
≤ 0
x
1
=0.000
x
2
=0.000
tak
6
0.045x
1
– x
2
≤ 0
0.4x
1
≤ 800
x
1
=2000.000
x
2
=90.000
nie
7
0.045x
1
– x
2
≤ 0
2.1x
1
+ 6x
2
≤ 600
x
1
=253.165
x
2
=11.392
tak
8
0.045x
1
– x
2
≤ 0
-x
1
≤ 0
x
1
=0.000
x
2
=0.000
tak
9
0.045x
1
– x
2
≤ 0
-x
2
≤ 0
x
1
=0.000
x
2
=0.000
tak
10
0.4x
1
≤ 800
2.1x
1
+ 6x
2
≤ 600
x
1
=2000.000
x
2
=-600.000
nie
11
0.4x
1
≤ 800
-x
1
≤ 0
- nie
12
0.4x
1
≤ 800
-x
2
≤ 0
x
1
=2000.000
x
2
=0.000
nie
13
2.1x
1
+ 6x
2
≤ 600
-x
1
≤ 0
x
1
=0.000
x
2
=100.000
nie
14
2.1x
1
+ 6x
2
≤ 600
-x
2
≤ 0
x
1
=285.714
x
2
=0.000
nie
15
-x
1
≤ 0
-x
2
≤ 0
x
1
=0.000
x
2
=0.000
tak
najlepsze rozwiązanie w punkcie x
1
=181.818; x
2
=36.364 (f(X)=236.363636 ).
Rozwiązanie 2
algorytm SIMPLEX rozwiązuje problemy w postaci
max
f(X)
przy
ograniczeniach
AX = b
X
≥ 0; b ≥ 0
(n – ilość zmiennych; m – ilość ograniczeń)
więc
max
{0.5x
1
+ 4x
2
}
x
1
, x
2
∈
Ω
-x
1
+ 5x
2
+ x
3
= 0
Ω
:
-x
2
+ 0.045x
1
+ x
4
= 0
0.4x
1
+ x
5
= 800
2.1x
1
+ 6x
2
+ x
6
= 600
x
j
≥ 0; j=1,2..n
n=6; m=4
m < n – problem posiada więcej niż jedno możliwe rozwiązanie
współczynniki przy zmiennych bazowych (x
j
; j=m+1..n) muszą być równe „1”. W przykładzie ze
względu na zerowe wartości parametru b w odpowiednich równaniach (1 i 2) możliwe było przemnożenie
przez „–1”. Jednak dla przykładowego ograniczenia
x
1
+ x
2
≥ 2
przekształcenie
x
1
+ x
2
– x
3
= 2 nie utworzy rozwiązania bazowego. W takiej
sytuacji tworzy się tzw. sztuczną bazę. Konieczne jest wprowadzenie dodatkowej zmiennej aby zachować
warunek b
≥ 0
Dla rozpatrywanego problemu
max {0.5x
1
+ 4x
2
}
x
1
, x
2
∈
Ω
x
1
- 5x
2
- x
3
+x
7
= 0
→
x
7
= -x
1
+ 5x
2
+ x
3
Ω
:
x
2
- 0.045x
1
- x
4
+ x
8
= 0
→
x
8
= -x
2
+ 0.045x
1
+ x
4
0.4x
1
+ x
5
= 800
2.1x
1
+ 6x
2
+ x
6
= 600
x
j
≥ 0; j=1,2..8
Wprowadza się dodatkową funkcje celu
T= x
7
+ x
8
=-0.955x
1
+ 4x
2
+ x
3
+ 4x
4
Tablica 1 Tablica SIMPEX dla rozpatrywanego przykładu bez wprowadzania sztucznej bazy
krok
obliczeń
nr ograniczenia
funkcja celu
b 1 2 3 4 5 6
b
i
/a
ij
; a
ij
>0
1 0
-1
5
1
0 0 0 =0/5
2 0
0.045
-1
0
1
0 0
-
3 800
0.4
0 0 0 1
0 -
4 600
2.1 6 0 0 0 1
=600/6
1
-f=- 0
-0.5
-4
0 0 0 0
-
1 0
-0.2
1
0.2 0 0 0
-
2 0
-0.155
0
0.2
1
0 0
-
3 800
0.4
0 0 0 1
0
=800/0.4
4 600
3.3
0 -1.2 0 0 1
=600/3.3
2
-f=- 0
-1.3
0 0.8 0 0 0
1
36.36
0
1
0.1273
0 0
0.0606
2 28.18
0 0
0.1436
1
0 0.0469
3 727.27 0
0 0.1454
0
1
-0.1212
rozwiązane
bazowe
4
181.81 1
0 -0.3636
0
0 0.3030
3
-f=- 236.36 0
0 0.3273
0
0 0.39393
Tablica 2 Tablica SIMPEX dla rozpatrywanego przykładu w wprowadzaniem sztucznej bazy
krok
obliczeń
nr ograniczenia
funkcja celu
b 1 2 3 4 5 6 7 8
b
i
/a
ij
; a
ij
>0
1 0
1
-5 -1 0 0 0 1
0
0/1
2
0 -0.045 1 0 -1 0 0 0 1
-
3 800
0.4
0
0
0
1
0 0 0 800/0.4
4
600
2.1 6 0 0 0 1
0 0 600/2.1
-f=-
0
-0.5
-4 0 0 0 0 0 0
1
T=
0
-0.955
4 1 1 0 0 0 0
1 0
1
-5 -1 0 0 0
0
-
2 0
0
0.775 -0.045
-1 0 0
1
0/0.775
3 800
0
2
0.4
0
1
0 0
800/2
4 600
0
16.5
2.1
0
0
1
0 600/16.5
-f=- 0
0
-6.5 -0.5 0 0 0
0
2
T=
0 0
-0.775 0.045
1 0 0
0
1 0
1
0 -1.29
-6.452
0 0
-
2 0
0
1
-0.058
-1.29 0 0
-
3 800
0
0
0.516
2.581
1
0 800/2.581
4 600
0
0
3.058 21.29
0
1
600/21.29
-f=-
0 0 0
-0.877 -8.387
0 0
3
T=
0 0 0 0 0 0 0
usuni
ęto z rozw. bazowego
x
7
usuni
ęto z rozw.
bazowego
x
8
II faza
1
181.81 1
0 -0.364
0 0 0.303
-
2
36.36
0
1
0.127 0 0 0.061
3
727.3 0 0 0.145 0 1
-0.121
4 28.18
0
0
0.144
1
0 0.047
4
-f=-
236.4 0 0 0.327 0 0 0.394
ujemny współczynnik
o największej wartości
bezwzględnej
element główny
najmniejszy
współczynnik
Skrypty Matlaba
wydruk skryptu wyliczającego wszystkie punkty przecięcia ograniczeń (przy 2 zmiennych decyzyjnych)
A=[-1 5
0.045
-1
0.4
0
2.1
6
-1
0
0
-1];
b=[0 0 800 600 0 0]';
x=[];
for i=1:size(A,1)-1
for j=i+1:size(A,1)
tmp=A([i j],:)\b([i j]);
x=[x tmp];
end
end
clear tmp i j
wydruk skryptu wykonującego pojedynczą iterację Gaussa-Jordana względem elementu głównego o
indeksie podanym w wywołaniu
%gj_1k.m - funkcja wykonujaca jedna iteracje gaussa-jordana wzgledem
%elementu glownego podanego w wywolaniu
%
% skladnia [MtrxOut]=gj_1k(Mtrx,IndxElGl)
%
%wejscie:
% Mtrx - macierz na ktorej wykonana bedzie iteracja Gaussa-Jordana
% IndxElGl - indeks elementu glownego w macierzy Mtrx
% IndxElGl(1) - wiersz, w ktorym znajduje sie element glowny
% IndxElGl(2) - kolumna, w ktorej znajduje sie element glowny
%
%wyjscie:
% MtrxOut - macierz Mtrx po wykonaniu jednej iteracji Gaussa-Jordanan na
% macierzy Mtrx wzgledem elementu o indeksie podanym w wywolaniu
%Autor : KD
%9 IV 2003 - utworzenie
function [MtrxOut]=gj_1k(Mtrx,IndxElGl);
% sprawdzenie poprawnosci danych wejsciowych
if nargin < 2, error('Za malo parametrow wejsciowych'); end
tmp=size(Mtrx);
ia=tmp(1);
ja=tmp(2);
if(ia*ja== length(Mtrx)|... %jesli zmienna z wywolania to wektor lub macierz ma wiecej niz 2 wymiary
length(tmp) > 2), error('Zmienna Mtrx musi byc macierza dwuwymiarowa'); end
if(length(IndxElGl(:)) <2) error('Zmienna IndxElGl musi byc wektorem'); end
if(IndxElGl(1) > ia | IndxElGl(1) < 1 |...
IndxElGl(2) > ja | IndxElGl(2) < 1 |...
sum(floor(IndxElGl)-IndxElGl) ~=0 ), error('nieprawidlowy indeks elementu glownego'); end
tmp=ones(ia,1)*(Mtrx(IndxElGl(1),:)/Mtrx(IndxElGl(1),IndxElGl(2)));
tmp= tmp .* (-Mtrx(:,IndxElGl(2))*ones(1,ja));
MtrxOut=Mtrx+tmp;
MtrxOut(IndxElGl(1),:)=(Mtrx(IndxElGl(1),:)/Mtrx(IndxElGl(1),IndxElGl(2)));
skrypt kreślący funkcje przestrzeń optymalizacji z ograniczeniami z rysunku 1
clc
clear
f=[.5 4];
A=[-1 5
0.045
-1
0.4
0
2.1
6
-1
0
0
-1];
b=[0 0 800 600 0 0]';
x1=0:100:2100;
x2=[-5:4:70]';
X1=ones(length(x2),1)*x1;
X2=x2*ones(1,length(x1));
F=f(1)*X1 + f(2)*X2;
mesh(X1,X2,F)
hold on;
colormap([0 0 1]);
x2_1=1/5*x1;
ogr=f(1)*x1 + f(2)*x2_1;
indx=find(x2_1 >= min(x2) & x2_1 <= max(x2));
plot3(x1(indx),x2_1(indx),ogr(indx),'k','LineWidth',1.5)
x2_1=0.045*x1;
ogr=f(1)*x1 + f(2)*x2_1;
indx=find(x2_1 >= min(x2) & x2_1 <= max(x2));
plot3(x1(indx),x2_1(indx),ogr(indx),'k','LineWidth',1.5)
x1_1=ones(length(x2),1)*800/0.4;
ogr=f(1)*x1_1 + f(2)*x2;
indx=find(x1_1 >= min(x1) & x1_1 <= max(x1));
plot3(x1_1(indx),x2(indx),ogr(indx),'k','LineWidth',1.5)
x2_1=( 600 - 2.1*x1)/6;
ogr=f(1)*x1 + f(2)*x2_1;
indx=find(x2_1 >= min(x2) & x2_1 <= max(x2));
plot3(x1(indx),x2_1(indx),ogr(indx),'k','LineWidth',1.5)
x1_1=zeros(length(x2),1);
ogr=f(1)*x1_1 + f(2)*x2;
plot3(x1_1,x2,ogr,'k','LineWidth',1.5)
x2_1=zeros(1,length(x1));
ogr=f(1)*x1 + f(2)*x2_1;
plot3(x1,x2_1,ogr,'k','LineWidth',1.5)
clear x1_1 x2_1 indx X1 X2
x=linprog(-f,A,b);
plot3(x(1), x(2), f*x,'ro','MarkerEdgeColor','r',...
'MarkerFaceColor','r',...
'MarkerSize',7);
x=[x [0 ; 0] (A([2 4],:)\b([2 4])), x];
ogr=f(1)*x(1,:) + f(2)*x(2,:);
plot3(x(1,:), x(2,:), ogr,'r','LineWidth',3);
axis([min(x1) max(x1) min(x2) max(x2) min(F(:)) max(F(:))]);
xlabel('x1')
ylabel('x2')
zlabel('f')
clear F x1 x2 ogr
x=x(:,1);