Rozdíly

Zde můžete vidět rozdíly mezi vybranou verzí a aktuální verzí dané stránky.

--- courses:a4b33opt:cviceni5 [2009/10/19 09:54]
r.polasek created
+++ courses:a4b33opt:cviceni5 [2025/01/03 18:28] (aktuální)
@@ Řádek 1: / Řádek 1: @@
-===== Cvičení 4 =====
+===== Cvičení 5 =====
+<code matlab>
+%figure;
+hold on;
+%primka
+a = rand(1);
+b = rand(1);
+xLine = (-5:5);
+yLine = a*xLine + b;
+plot(xLine, yLine, 'y--')
+%body v normalnim rozdeleni
+epsilon = 0.05;
+v = randn(2, 100);
+sx = v(1, 1:length(v));
+x = sx(:);
+y = v(2, 1:length(v));
+y = a * x + b+ y(:);
+plot(x, y, 'r+')
+%robustni prokladani primarni ulohy
+heigth = length(x);
+y1 = y+epsilon;
+y2 = -y+epsilon;
+y3 = 0;
+A = [ x ones(heigth,1) -1*eye(heigth); -1*x -1*ones(heigth,1) -1*eye(heigth); zeros(heigth,1) zeros(heigth,1) -1*eye(heigth) ];
+b = [y1;y2;zeros(heigth,1)];
+f = [0; 0; ones(heigth,1)];
+v = linprog(f, A, b );
+PrimaryOptimum = f'*v
+plot(x,y,'r+');
+lineX = [-5:0.1:5];
+lineX = lineX(:);
+lineY = v(1)*lineX + v(2);
+plot(lineX,lineY, 'b');
+plot(lineX,lineY+epsilon, 'g--');
+plot(lineX,lineY-epsilon, 'g--');
+hold off;
+%robustni prokladani dualni ulohy
+fD = -b';
+AD = A';
+bD = f';
+AeqD = (A'); %tady skromnej trik...nema to bejt omezeny na nerovnost, ale na rovnost....
+AeqD = AeqD(1:2,:); %protoze jsem linej psat znova ty nerovnice (aD, bD), tak je tam necham, jak jsou.
+beqD = f'; %Protoze je rovnost tvrdsi omezeni, tak se jejim pridanim nic neporusi a docilim
+beqD = beqD(1:2); %stejnych vlasnosti, jako kdyby aD, bD ony dve rovnice nebyly a byly jen v AeqD a beqD
+v2 = linprog(fD, AD, bD,AeqD,beqD,[],zeros(length(b), 1));
+SecondaryOptimum = -fD*v2
+%vykresleni divergenci a odpovidajicich dualnich promennych
+figure
+hold on
+xAxis = [1:1:length(v2)];
+plot(xAxis, v2, 'r+')
+hold off;
+</code>
+==== Bonus 1 ====
+ --- //[[[email protected]|Marek Sacha]] 2009/10/29 17:16//
+<code matlab>
+clear
+clc
+load data1.mat
+% Construct inequalities matrix
+%       | x   1  -1 |
+% A =   |-x  -1  -1 |
+%
+A = [ x ones(50,1) -1*ones(50,1); -1*x -1*ones(50,1) -1*ones(50,1) ];
+% Construct right side
+%     |  y |
+% b = | -y |
+b = [ y; -y ];
+% Construct function to be minimized
+% f = 0x + 0y + lambda
+f = [ 0; 0; 1];
+% Run linprog
+o = linprog(f,A,b);
+% o = [ a b lambda ];
+% Modify linprog parameters for secondary problem
+% Run linprog
+f2 = b';
+A2 = -1*A';
+b2 = f';
+[ awh a2w ] = size(A2);
+[ f2h f2w ] = size(f2);
+lb = zeros(f2h,f2w);
+% Dont forget to add equality conditions (for first two rows of A2,b2)
+% - a, b are from 'R' not R+
+o2 = linprog(f2,A2,b2,A2(1:2,[1:a2w]),b2(1:2),lb);
+% Verify results
+result_prim = (-f)'*o
+result_sec = b'*o2
+hold on;
+% Print separate points
+plot(x,y,'y.');
+% Approximation line
+plot(x,x*o(1)+o(2),'r');
+% Top border
+plot(x,x*o(1)+o(2)+o(3),'g');
+% Bottom border
+plot(x,x*o(1)+o(2)-o(3),'g');
+hold off;
+% New Plot
+figure;
+% Visualize dual
+plot([x ; x],[y ; -y].*o2,'ro');
+</code>
+==== Bonus 2 ====
+ --- //[[[email protected]|Marek Sacha]] 2009/10/29 17:16//
+<code matlab>
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Author: Marek Sacha <[email protected]> %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+clear all
+clc
+numberOfPoints = [];
+primal = [];
+dual = [];
+options = optimset('LargeScale','off','Simplex','on','MaxIter',1e5);
+for j=1:4
+    t = j*100;
+    realloc = [ numberOfPoints t ];
+    numberOfPoints = realloc;
+    % initialize random matrix
+    r = rand(t,2);
+    x = r(:,1);
+    y = r(:,2);
+    % Initialize e
+    e = 0.05;
+    [ xh xw ] = size(x);
+    % Construct inequalities matrix
+    %         a    b   l1  l2 ... ln
+    %        ++++++++++++++++++++++++
+    %       | x1   1  -1   0  ...  0 |
+    %       | x2   1   0  -1  ...  0 |
+    %       ...
+    %       | xn   1   0   0  ... -1 |
+    % A =   |-x1  -1  -1   0  ...  0 |
+    %       |-x2  -1   0  -1  ...  0 |
+    %       ...
+    %       | 0    0   0   0  ... -1 |
+    %        ++++++++++++++++++++++++
+    %
+    A = [ x ones(xh,1) -1*eye(xh); -1*x -1*ones(xh,1) -1*eye(xh); zeros(xh,1) zeros(xh,1) -1*eye(xh) ];
+    % Construct right side
+    %     | e + y |
+    % b = | e - y |
+    %     | 0     |
+    b = [ (e*ones(xh,1))+y; (e*ones(xh,1))-y; zeros(xh,1) ];
+    % Construct function to be minimized
+    % f = 0a + 0b + lambda1 + lambda2 + lambda3 + ... + lambdan
+    f = [ 0 ; 0; ones(xh,1)];
+    % Run linprog
+    tic;
+    o = linprog(f,A,b,[],[],[],[],[],options);
+    time = toc;
+    realloc = [ primal time ];
+    primal = realloc;
+    % Modify linprog parameters for secondary problem
+    f2 = b';
+    A2 = -1*A';
+    b2 = f';
+    [ awh a2w ] = size(A2);
+    [ f2h f2w ] = size(f2);
+    lb = zeros(f2h,f2w);
+    % Dont forget to add equality conditions (for first two rows of A2,b2)
+    % - a, b are from 'R' not R+
+    tic;
+    o2 = linprog(f2,A2,b2,A2(1:2,[1:a2w]),b2(1:2),lb,[],[],options);
+    time = toc;
+    realloc = [ dual time ];
+    dual = realloc;
+    % Verify results
+    result_prim = (-f)'*o;
+    result_sec = b'*o2;
+end
+hold on;
+plot(numberOfPoints,primal,'g-');
+plot(numberOfPoints,dual,'r-');
+hold off;
+</code>
 ~~DISCUSSION~~

courses/a4b33opt/cviceni5.1255938842.txt.gz · Poslední úprava: 2025/01/03 18:23 (upraveno mimo DokuWiki)

Nahoru