%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;

— Marek Sacha 2009/10/29 17:16

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

— Marek Sacha 2009/10/29 17:16

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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;