Ć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);