Метод динамического программирования

2.1 Теоретическая часть

Математически задачу выбора набора параметров из заданной их совокупности можно сформулировать следующим образом.

Пусть работоспособность объекта контроля характеризуется совокупностью n взаимосвязанных параметров. Образующих множество S={x₁, x₂,…,x_n}. Проверка всех параметров из S влечет контроль всех N элементов системы и дает однозначный ответ: объект исправен, если все N элементов исправны, или неисправен. Если по крайней мере один из элементов отказал. Для x_i определено подмножество R(x_i) элементов, проверяемых при контроле i-ого параметра. Причем предполагается, что эти подмножества могут пересекаться, т.е. i, j: R(x_i) ∩ R(x_j).

Пусть Ω – некоторый набор параметров из множества S, т.е. . Тогда и , где Ω-набор контролируемых параметров, а - неконтролируемый набор. Значение из S можно представить булевым вектором, причем:

1, если xi

0. если xi

xi=

Задача выбора параметров в этом случае формулируется двояко:

1) Найти набор, для которого

Р(Ω)=max (1.11)

при

2) Найти набор Ω. Для которого:

(1.12)

при

где - апостериорная вероятность работоспособного состояния объекта контроля при положительном исходе контроля выбранных параметров;

c{x_i} - затраты на контроль i-ого параметра;

- требуемая достоверность контроля;

С – ограничение на общую стоимость контроля.

Значение зависит от принятых допущений и может быть найдено по формуле

Байеса. Так, если полагать в изделии наличие лишь одного отказа, то

(1.13)

где - априорная вероятность безотказной работы объекта:

- нормированная вероятность отказа системы из-за отказа i-ого элемента:

(1.14)

- априорная вероятность отказа i-ого элемента.

Тогда вероятность того, что отказ будет обнаружен при проверки i-ого параметра можно вычислить по формуле:

При возможности наличия в ОК произвольного числа отказов

(1.15)

Для решения задач (1.11) и (1.12) можно использовать простой перебор вариантов, однако возникающие при этом вычислительные трудности не позволяют сделать этого даже для простых систем ( при n>10). В связи с этим комплектование набора трактуется как многошаговый процесс, состоящий из последовательного выбора отдельных параметров. В соответствии с общим принципом оптимальности, разобьем весь имеющийся ресурс стоимости С на С отрезков единичной длины. (В практических случаях заданные положительные величины c(x_i) и С можно считать всегда целыми. Если это не так, то необходимо перейти к более мелким стоимостным единицам в зависимости от разрядности дробной части.) Рассмотрим наряду с интересующей нас исходной задачей множество аналогичных задач

, аналог выражения (1.1),

где через X_С обозначено множество неотрицательных целочисленных векторов Ω, отвечающих наборам, в которых общая стоимость проверки параметров не превосходит величины С.

Пусть

тогда при всех соответствующие множества X_Y состоят. Очевидно, из одного нулевого элемента и f(С)=0 для всех таких С. Для ресурса согласно общей схеме динамического программирования справедливы следующие рекуррентные соотношения:

(1.16)

где

R( )-множество тех i, для которых , начиная с номера уравнение (1.16) решается для всех i=1,n:

- сумма вероятностей элементов i-го параметра, которые пересекаются с элементами подмножества , образованного на шаге

Если

Æ, то

И

(1.17)

Что приводит к условию задачи с пересекающимися параметрами.

Для решения рассматриваемой задачи рассмотрим простой численный метод, не требующий предварительного определения всех допустимых наборов и основанный на рекуррентных соотношениях (1.15). Для всех целых по формуле (1.16) вычисляют величины и при этом фиксируются индексы , на которых достигаются максимумы в (1.16). Искомый вектор формируется последовательно включение в набор параметра и подмножества , зафиксированного на шаге При этом, если , то на данном шаге этот параметр исключается из рассмотрения, так как каждый параметр может включаться в набор не более одного раза. Если на н6екотором v-м шаге окажется, что то вы качестве принимается я подмножество и фиксируется параметр , причем за принимается значение . Заметим, что если в (1.11) принять более жесткое ограничение, а именно

то последнее недопустимо, так как в этом случае может быть меньше из-за того, что он достигает на и другом подмножестве параметров.

2.2 Практическая часть

Дано:

Нормированные вероятности отказов для элементов Q_i, затраты на контроль i-ого параметра C(x_i).

№
Qi	0,04	0,15	0,13	0,08	0,02	0,09	0,06	0,08	0,03	0,04
C(xi)

Найти: Выбрать такие параметры, чтобы C≤18 при Q=Q_max.

Для удобства расчетов проранжируем таблицу следующим образом:

№
Qi	0,15	0,04	0,03	0,08	0,06	0,09	0,08	0,13	0,04	0,02
C(xi)

Решение:

Весь имеющийся ресурс разбиваем на отрезки единичной длины. Таким образом, .

При соответствующее множество состоит из одного нулевого элемента и целевая функция .

При : в множество можем включить один из трёх подходящих параметров по стоимости затрат, но выбираем тот параметр, который обладает наибольшей вероятностью, а именно , вероятность которого равна . Вычислили значение целевой функции и зафиксировали индекс параметра .

При : на этом шаге для включения в множество можем выбрать параметр со стоимостью, равной единицы, тогда множество будет состоять из двух параметров, или же можем выбрать со стоимостью, равной два. Соответственно выбираем c вероятностями . Целевая функция

и фиксируем .

Аналогично проводим расчёты на последующих шагах. Результат пошаговой работы алгоритма представлен в виде таблицы 2.

Таблица 2.

Выделяемое количество затратных единиц, Ck	Значение целевой функции, f(Ck)	Выбор на данном шаге, i*Ck	Выбранный массив, Ω*
		-
	0,15
	0,19		2;10
	0,23		2;4
	0,27		2;4;10
	0,3		2;4;10;9
	0,33		2;4;10;7
	0,36		2;4;10;7;9
	0,4		2;4;10;3
	0,43		2;4;10;3;9
	0,44		2;4;10;3;7
	0,47		2;4;10;3;7;9
	0,52		2;4;10;3;9;6
	0,55		2;4;10;3;6;7
	0,58		2;4;10;3;6;7;9
	0,6		2;4;10;3;6;9;8
	0,63		2;4;10;3;6;8;7
	0,66		2;4;10;3;6;8;7;9
	0,66	-	2;4;10;3;6;8;7;9

Таким образом, решением является набор проверок Ω* = {X₂;X₃;X₄;X₆;X₇;X₈;X₉;X₁₀}.

2.3 Программные расчеты и сравнение результатов

Наряду с ручным расчётом, решение задачи реализовано с помощью программного алгоритма, написанного на языке программирования Delphi 7.0.

Листинг программы представлен в приложении 2.

Результаты программного расчёта сохраняются в текстовый файл kr2.txt и представлены на рисунке 3.

Рис.3. Результат работы программы для метода динамического программирования

Выводы

При решении поставленной задачи оба метода дали одинаковый результат, что показывает правильность выполнения и, соответственно, понимания курсовой работы. Таким образом, необходимо использовать следующие параметры контроля: .

Метод ветвей и границ является вариацией метода полного перебора с той разницей, что мы сразу исключаем заведомо неоптимальные решения. Главный недостаток алгоритма метода ветвей и границ заключается в необходимости полностью решать задачи линейного программирования, ассоциированные с каждой из вершин многогранника допустимых решений. Для задач большой размерности это требует значительных затрат времени.

Динамическое программирование – один из наиболее мощных методов оптимизации. Этот метод исключительно привлекателен благодаря простоте и ясности своего основного принципа – принципа оптимальности.

Этот метод отличается от линейного программирования, исходная постановка основных задач которых имеет статический характер. Главным недостатком метода является – его сложность возрастает с увеличением размерности задачи.

Список используемой литературы

1. Пискунов В.А. Курс лекций и лабораторных работ по дисциплине “Контроль и диагностика систем ЛА”.

2. Пискунов В.А. Учебное пособие “Комбинаторные методы дискретной оптимизации в прикладных задачах организации контроля систем ЛА”.

Приложение 1

Листинг программы

program vetvi;

const

maxmatrix=1000;

maxsize=200;

type linear=array[1..10000] of integer;

label skip;

var matrix:array[1..1000] of pointer;

n:integer;

sizeofm:word;

q,w,e,r:integer;

start_m:integer;

sm:^linear;

bestx,besty:integer;

bz:integer;

ochered:array[1..1000] of

record

id:integer;

ocenka:integer;

end;

nochered:integer;

workm,workc:integer;

leftm,rightm:integer;

first,last:integer;

best:integer;

bestmatr:array[1..maxsize] of integer;

bestmatr1:array[1..maxsize] of integer;

curr:integer;

procedure swapo(a,b:integer);

begin

ochered[1000]:=ochered[a];

ochered[a]:=ochered[b];

ochered[b]:=ochered[1000];

end;

procedure addochered(id,ocenka:integer);

var

curr:integer;

begin

inc(nochered);

ochered[nochered].id:=id;

ochered[nochered].ocenka:=ocenka;

{Uravnoveshivanie ocheredi}

curr:=nochered;

while true do

begin

if curr=1 then break;

if ochered[curr].ocenka< ochered[curr div 2].ocenka

then

begin

swapo(curr,curr div 2);

curr:=curr div 2;

end

else break;

end;

end;

procedure getochered(var id,ocenka:integer);

var

curr:integer;

begin

id:=ochered[1].id;

ocenka:=ochered[1].ocenka;

ochered[1]:=ochered[nochered];

dec(nochered);

curr:=1;

while true do

begin

if (curr*2+1> nochered) then break;

if (ochered[curr*2].ocenka< ochered[curr].ocenka) or

(ochered[curr*2+1].ocenka<ochered[curr].ocenka) then

begin

if ochered[curr*2].ocenka> ochered[curr*2+1].ocenka

then

begin

swapo(curr*2+1,curr);

curr:=curr*2+1;

end

else

begin

swapo(curr*2,curr);

curr:=curr*2;

end;

end else break;

end;

end;

function getid:integer;

var q:integer;

qw:^linear;

begin

if memavail<10000 then

begin

q:=ochered[nochered].id;

{ exit;}

end else

begin

for q:=1 to maxmatrix do

if matrix[q]=nil then break;

getmem(matrix[q],sizeofm);

end;

qw:=matrix[q];

fillchar(qw^,sizeofm,0);

getid:=q;

end;

procedure freeid(id:integer);

begin

freemem(matrix[id],sizeofm);

matrix[id]:=nil;

end;

function i(x,y:integer):integer;

begin

i:=(y-1)*n+x+1;

end;

function simplize(id:integer):integer;

var

q,w:integer;

t:^linear;

add:integer;

min:integer;

begin

t:=matrix[id];

add:=0;

for q:=1 to n do

begin

min:=maxint;

for w:=1 to n do

if t^[i(w,q)]< >-1 then

if min> t^[i(w,q)] then min:=t^[i(w,q)];

if min<>0 then

for w:=1 to n do

if t^[i(w,q)]< >-1 then

dec(t^[i(w,q)],min);

if min>32000 then min:=0;

inc(add,min);

end;

for q:=1 to n do

begin

min:=maxint;

for w:=1 to n do

if t^[i(q,w)]< >-1 then

if min> t^[i(q,w)] then min:=t^[i(q,w)];

if min<>0 then

for w:=1 to n do

if t^[i(q,w)]< >-1 then

dec(t^[i(q,w)],min);

if min>32000 then min:=0;

inc(add,min);

end;

simplize:=add;

end;

function bestziro(id:integer):integer;

var

t:^linear;

q,w,e,x,y:integer;

min1,min2:integer;

l1,l2:array[1..maxsize] of integer;

begin

t:=matrix[id];

fillchar(l1,sizeof(l1),0);

fillchar(l2,sizeof(l2),0);

for q:=1 to n do

begin

min1:=maxint;min2:=maxint;

for w:=1 to n do

if t^[i(w,q)]< >-1 then

begin

if min2> t^[i(w,q)] then min2:=t^[i(w,q)];

if min1> min2 then

begin

e:=min1;

min1:=min2;

min2:=e;

end;

end;

if min1<>0 then min2:=0;

if min2>32000 then min2:=0;

l2[q]:=min2;

end;

for q:=1 to n do

begin

min1:=maxint;min2:=maxint;

for w:=1 to n do

if t^[i(q,w)]< >-1 then

begin

if min2> t^[i(q,w)] then min2:=t^[i(q,w)];

if min1> min2 then

begin

e:=min1;

min1:=min2;

min2:=e;

end;

end;

if min1<>0 then min2:=0;

if min2>32000 then min2:=0;

l1[q]:=min2;

end;

bz:=-32000;

bestx:=0;besty:=0;

for y:=n downto 1 do

for x:=1 to n do

if (t^[i(x,y)]=0) then

if l1[x]+l2[y]> bz then

begin

bestx:=x;

besty:=y;

bz:=l1[x]+l2[y];

end;

bestziro:=bz;

end;

begin

assign(input,'input.txt');

assign(output,'vetvi.out');

reset(input);

rewrite(output);

nochered:=0;

read(n);

sizeofm:=n*(n+2)*2+2;

start_m:=getid;

sm:=matrix[start_m];

for q:=1 to n do

for w:=1 to n do

read(sm^[i(w,q)]);

addochered(start_m,0);

{ ;

bestziro(start_m);}

{Sobstvenno reshenie}

best:=maxint;

while true do

begin

if nochered=0 then break;

getochered(workm,workc);

{process MATRIX}

inc(workc,simplize(workm));

if workc> best then goto skip;

sm:=matrix[workm];

if sm^[1]=n-1 then

begin

best:=workc;

for q:=1 to n do

begin

bestmatr [q]:=sm^[i(q,n+2)];

bestmatr1[q]:=sm^[i(q,n+1)];

end;

goto skip;

end;

q:=bestziro(workm);

if q=-32000 then goto skip;

{Pravaia vetka}

if(bestx=0) or (besty=0) then goto skip;

rightm:=getid;

move(matrix[workm]^,matrix[rightm]^,sizeofm);

sm:=matrix[rightm];

sm^[i(bestx,besty)]:=-1;

addochered(rightm,workc+q);

{Levaia vetka}

leftm:=getid;

move(matrix[workm]^,matrix[leftm]^,sizeofm);

sm:=matrix[leftm];

{Dobavliaetsia rebro iz bestx v besty}

inc(sm^[1]);

sm^[i(bestx,n+2)]:=besty;

sm^[i(besty,n+1)]:=bestx;

first:=bestx;last:=besty;

if sm^[1]< >n-1 then

begin

while true do

begin

if sm^[i(last,n+2)]=0 then break;

last:=sm^[i(last,n+2)];

end;

while true do

begin

if sm^[i(first,n+1)]=0 then break;

first:=sm^[i(first,n+1)];

end;

sm^[i(last,first)]:=-1;

sm^[i(first,last)]:=-1;

sm^[i(besty,bestx)]:=-1;

end;

for w:=1 to n do

begin

sm^[i(w,besty)]:=-1;

sm^[i(bestx,w)]:=-1;

end;

addochered(leftm,workc);

skip:

{Free Matrix}

freeid(workm);

end;

{ freeid(start_m);}

if best=maxint then

begin

writeln('Путь не существует');

end else

begin

writeln('Длина пути:',best);

for q:=1 to n do

if bestmatr[q]=0 then break;

e:=q;

for curr:=1 to n do

if bestmatr[curr]=q then break;

while true do

begin

write(curr,' ');

curr:=bestmatr1[curr];

if curr=0 then

begin

writeln(e);

break;

end;

end;

end;

close(input);

close(output);

end.

Приложение 2

Листинг программы

unit Unit1;

interface

uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

Dialogs, ToolWin, ComCtrls, mdCONTROLS, Grids, StdCtrls, ExtCtrls, Unit2,

Buttons;

type

TForm1 = class(TForm)

sgH: TStringGrid;

sgP: TStringGrid;

sgC: TStringGrid;

sgQ: TStringGrid;

lbC: TLabeledEdit;

BitBtn1: TBitBtn;

Label1: TLabel;

sgW: TStringGrid;

Label2: TLabel;

procedure FormCreate(Sender: TObject);

procedure BitBtn1Click(Sender: TObject);

procedure sgExit(Sender: TObject);

private

{ Private declarations }

public

H: TH;

P: TP;

C: TC;

W: TW;

end;

var

Form1: TForm1;

implementation

{$R *.dfm}

procedure TForm1.FormCreate(Sender: TObject);

var i,j: integer;

x: Byte;

f: TextFile;

begin

AssignFile(f, 'data.txt');

Reset(f);

sgW.Cells[0,0] := 'W';

// Ввод исходной матрицы

readln(f);

for j:=1 to 10 do

begin

sgH.Cells[0,j] := IntToStr(j);

sgW.Cells[0,j] := IntToStr(j);

for i:=1 to 10 do

begin

sgH.Cells[i,0] := IntToStr(i);

read(f, x);

sgH.Cells[i,j] := IntToStr(x);

if x = 1 then

H[i-1,j-1] := true

else

H[i-1,j-1] := false;

end;

readln(f);

end;

// Ввод вероятностей

readln(f);

readln(f);

sgP.Cells[0,0] := 'P';

for i:=1 to 10 do

begin

read(f, P[i-1]);

sgP.Cells[i,0] := FloatToStr(P[i-1]);

end;

readln(f);

// Ввод стоимостей

readln(f);

readln(f);

sgC.Cells[0,0] := 'C';

for j:=1 to 10 do

begin

read(f, C[j-1]);

sgC.Cells[0,j] := IntToStr(C[j-1]);

end;

CloseFile(f);

// Ввод вероятностей обнаружения отказа

sgQ.Cells[0,0] := 'Q';

for j:=1 to 10 do

sgQ.Cells[0,j] := FloatToStr(Q(j-1,H,P));

lbC.Text := '1';

end;

procedure TForm1.BitBtn1Click(Sender: TObject);

var i: integer;

begin

label1.Caption := FloatToStr(maxf(1, StrToInt(lbC.Text), H,P,C, W));

for i:=1 to 10 do

begin

sgW.Cells[2,i] := IntToStr(W[i-1].N);

if W[i-1].E then

sgW.Cells[1,i] := '1'

else

sgW.Cells[1,i] := '0';

end;

end;

procedure TForm1.sgExit(Sender: TObject);

var i,j: integer;

begin

for j:=1 to 10 do

for i:=1 to 10 do

if sgH.Cells[i,j] = '1' then

H[i-1,j-1] := true

else

H[i-1,j-1] := false;

for i:=1 to 10 do

P[i-1] := StrToFloat(sgP.Cells[i,0]);

for j:=1 to 10 do

C[j-1] := StrToInt(sgC.Cells[0,j]);

// Ввод вероятностей обнаружения отказа

for j:=1 to 10 do

sgQ.Cells[0,j] := FloatToStr(Q(j-1,H,P));

end;

end.

unit Unit2;

interface

type

TH = array [0..9, 0..9] of boolean;

TP = array [0..9] of extended;

TC = array [0..9] of integer;

TDateW = record

E: boolean;

N: integer;

end;

TW = array [0..9] of TDateW;

function Q(j: integer; H: TH; P: TP): extended;

function maxf(n, Yk: integer; H: TH; P: TP; C: TC; var W: TW): extended;

implementation

function Q(j: integer; H: TH; P: TP): extended;

var i: integer;

begin

Result := 0;

for i:=0 to 9 do

if H[i,j] then

Result := Result + P[i];

end;

function G(j: integer; H: TH; P: TP; W: TW): extended;

var i,k: integer;

begin

Result := 0;

for i:=0 to 9 do

if H[i,j] then

for k:=0 to 9 do

if W[k].E and H[i,k] then

begin

Result := Result + P[i];

Break;

end;

end;

function f(n, Yk, j: integer; H: TH; P: TP; C: TC; var W: TW): extended;

begin

Result := Q(j,H,P) + maxf(n+1, Yk - C[j], H,P,C, W) - G(j,H,P,W);

end;

function maxf(n, Yk: integer; H: TH; P: TP; C: TC; var W: TW): extended;

var j,i: integer;

ft: extended;

Wt: TW;

begin

Result := 0;

for i:=0 to 9 do

begin

W[i].E := false;

W[i].N := 0;

end;

for j:=0 to 9 do

if C[j] <= Yk then

begin

for i:=0 to 9 do

begin

Wt[i].E := false;

Wt[i].N := 0;

end;

ft := f(n, Yk, j, H,P,C, Wt);

if Result < ft then

begin

Result := ft;

W := Wt;

W[j].E := true;

W[j].N := n;

end;

end;

end;

end.