[37] | 1 | clc |
---|
| 2 | echo on |
---|
| 3 | |
---|
| 4 | % This examle will show how to define LP and QP problems |
---|
| 5 | % |
---|
| 6 | % We will use YALMIP to formulate couple of regression problems |
---|
| 7 | % (linear L_1, L_2 and L_inf regression) |
---|
| 8 | % |
---|
| 9 | pause % Strike any key to continue. % Strike any key to continue. |
---|
| 10 | |
---|
| 11 | % To begin with, we generate data according to a noisy linear regression |
---|
| 12 | % |
---|
| 13 | % y(t) = x'(t)*a+e(t), t=0,0.01,...,2*pi |
---|
| 14 | |
---|
| 15 | a = [1 2 3 4 5 6]; |
---|
| 16 | t = (0:0.02:2*pi)'; |
---|
| 17 | x = [sin(t) sin(2*t) sin(3*t) sin(4*t) sin(5*t) sin(6*t)]; |
---|
| 18 | y = x*a'+(-4+8*rand(length(x),1)); |
---|
| 19 | pause % Strike any key to continue. |
---|
| 20 | plot(y) |
---|
| 21 | pause % Strike any key to continue. |
---|
| 22 | |
---|
| 23 | % We now want to estimate the parameter a, given y, and we will |
---|
| 24 | % try 3 different approaches; L_1, L_2 and L_inf. |
---|
| 25 | % |
---|
| 26 | % We start by defining the decision variable |
---|
| 27 | |
---|
| 28 | a_hat = sdpvar(1,6); |
---|
| 29 | pause % Strike any key to continue. |
---|
| 30 | |
---|
| 31 | % By using a_hat and the regressor x, we can define the residuals |
---|
| 32 | |
---|
| 33 | residuals = y-x*a_hat'; |
---|
| 34 | pause % Strike any key to continue. |
---|
| 35 | |
---|
| 36 | % To solve the L_1 problem, we define a variable that will work |
---|
| 37 | % as a bound on |y-x'*a| |
---|
| 38 | bound = sdpvar(length(residuals),1); |
---|
| 39 | pause % Strike any key to continue. |
---|
| 40 | |
---|
| 41 | % Express that bound is absolute value of the residuals (usinhg a double-sided constraint) |
---|
| 42 | F = set(-bound < residuals < bound); |
---|
| 43 | pause % Strike any key to continue. |
---|
| 44 | |
---|
| 45 | % Minimize sum of bound |
---|
| 46 | % NOTE, This takes time if you only have linprog installed |
---|
| 47 | solvesdp(F,sum(bound)); |
---|
| 48 | % Optimal L_1 estimate |
---|
| 49 | a_L1 = double(a_hat) |
---|
| 50 | pause % Strike any key to continue. |
---|
| 51 | |
---|
| 52 | |
---|
| 53 | clc |
---|
| 54 | |
---|
| 55 | % The L_2 solution is easily obtained from a QP |
---|
| 56 | % (note the nonlinear objetive function. More about |
---|
| 57 | % this in demo #11) |
---|
| 58 | solvesdp([],residuals'*residuals); |
---|
| 59 | a_L2 = double(a_hat) |
---|
| 60 | pause % Strike any key to continue. |
---|
| 61 | |
---|
| 62 | % Finally, let us minimize the L_inf norm. This corresponds to |
---|
| 63 | % minimizing the largest (absolute value) residual |
---|
| 64 | % We introduce a scalar bound on the largest value |
---|
| 65 | bound = sdpvar(1,1); |
---|
| 66 | F = set(-bound < residuals < bound); |
---|
| 67 | pause % Strike any key to continue. |
---|
| 68 | |
---|
| 69 | % Solve by minimizing the bound |
---|
| 70 | solvesdp(F,bound); |
---|
| 71 | a_Linf = double(a_hat) |
---|
| 72 | pause % Strike any key to continue. |
---|
| 73 | |
---|
| 74 | plot([y x*a' x*a_L1' x*a_L2' x*a_Linf']) |
---|
| 75 | legend('Measured','True','L1','L2','Linf') |
---|
| 76 | |
---|
| 77 | echo off |
---|
| 78 | |
---|
| 79 | |
---|
| 80 | |
---|
| 81 | |
---|
| 82 | |
---|
| 83 | |
---|
| 84 | |
---|
| 85 | |
---|
| 86 | |
---|
| 87 | |
---|