最优化理论作业.docx

Homework for CVX optimization
Name: **** Student ID : ********
For do this homework,I download the latest version Cvx Matlab toolbox, and get an academic license from CVX Research, Inc. All the programs below are completed with this toolbox, and can running under Matlab R2014a.
1 Example (Page 296)
We take a matrix A∈R100×30 and vector b∈R100 (chosen at random, but the results are typical), and compute the ℓ1-norm and ℓ2-norm approximate solutions of Ax≈b, as well as the penalty function approximations with a deadzone-linear penalty (with a = 0.5) and log barrier penalty (with a = 1). Figure 6.2 shows the four associated penalty functions,and the amplitude distributions of the optimal residuals for these four penalty approximations.
Solution:
I choose the matrix A for all 4 problems, and vector b with the former 3 problems. When I try to solve the last problem with the same b, the problem becomes to “Infeasible”. The reason is that most value of Ax-b are larger than 1.After trying many times, finally, I generate vector b using the randn () function, and divide by 1.5.
1) For ℓ1-norm, the object function is Ax-b , and using the following code to solve the problem
cvx_begin
variable x(n);
minimize( norm(A*x-b,1) );
cvx_end
The running result is as following:
Calling SDPT3 4.0: 230 variables, 100 equality constraints
------------------------------------------------------------
num. of constraints = 100
dim. of socp var = 200, num. of socp blk = 100
dim. of free var = 30 *** convert ublk to lblk
*******************************************************************
SDPT3: Infeasible path-following algorithms
*******************************************************************
version predcorr gam expon scale_data
NT 1 0.000 1 0
it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime
-------------------------------------------------------------------
0|0.000|0.000|9.2e-01|5.9e+01|3.0e+05| 2.497247e+02 0.000000e+00| 0:0:00| chol 1 1
1|1.000|0.882|9.9e-06|7.0e+00|1.9e+04| 4.908632e+03 5.514004e+01| 0:0:00| chol 1 1
2|1.000|0.990|3.9e-06|9.1e-02|1.5e+03| 1.491205e+03 2.675876e+01| 0:0:00| chol 1 1
3|0.945|1.000|2.6e-06|3.0e-03|8.2e+01| 1.133779e+02 3.184938e+01| 0:0:00| chol 1 1
4|0.817|0.203|3.2e-05|2.5e-03|3.5e+01| 7.128134e+01 3.697355e+01| 0:0:00| chol 1 1
5|0.898|0.243|7.9e-06|1.9e-03|2.2e+01| 6.403791e+01 4.243380e+01| 0:0:00| chol 1 1
6|0.860|0.441|9.0e-06|1.0e-03|1.2e+01| 6.165978e+01 4.996343e+01| 0:0:00| chol 1 1
7|0.905|0.518|3.3e-06|5.1e-04|5.3e+00| 6.025036e+01 5.499107e+01| 0:0:00| chol 1 1
8|0.929|0.526|4.8e-07|2.4e-04|2.4e+00| 5.984393e+01 5.745881e+01| 0:0:00| chol 1 1
9|1.000|0.508|5.9e-08|1.2e-04|1.1e+00| 5.972690e+01 5.859379e+01| 0:0:00| chol 1 1
10|1.000|0.172|2.4e-08|1.1e-04|9.7e-01| 5.973848e+01 5.877962e+01| 0:0:00| chol 1 1
11|1.000|0.468|9.5e-09|6.1e-05|5.3e-01| 5.972730e+01 5.920081e+01| 0:0:00| chol 1 1
12|0.982|0.754|3.4e-09|1.5e-05|1.3e-01| 5.969978e+01 5.957286e+01| 0:0:00| chol 1 1
13|0.885|0.288|8.6e-10|1.1e-05|9.1e-02| 5.969795e+01 5.960793e+01| 0:0:00| chol 1 1
14|1.000|0.269|3.8e-10|7.8e-06|6.7e-02| 5.969739e+01 5.963124e+01| 0:0:00| chol 1 1
15|0.859|0.771|1.8e-10|1.8e-06|1.5e-02| 5.969590e+01 5.968070e+01| 0:0:00| chol 1 1
16|0.983|0.966|2.1e-11|2.1e-05|6.8e-04| 5.969550e+01 5.969498e+01| 0:0:00| chol 1 1
17|0.989|0.989|2.4e-13|9.2e-07|1.1e-05| 5.969549e+01 5.969549e+01| 0:0:00| chol 1 1
18|0.833|0.945|4.1e-14|1.5e-08|4.6e-07| 5.969549e+01 5.969549e+01| 0:0:00|
stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
number of iterations = 18
primal objective value = 5.96954941e+01
dual objective value = 5.96954938e+01
gap := trace(XZ) = 4.62e-07
relative gap = 3.84e-09
actual relative gap = 3.25e-09
rel. primal infeas = 4.14e-14
rel. dual infeas = 1.46e-08
norm(X), norm(y), norm(Z) = 1.3e+01, 9.0e+00, 1.3e+01
norm(A), norm(b), norm(C) = 7.8e+01, 1.2e+01, 1.1e+01
Total CPU time (secs) = 0.29
CPU time per iteration = 0.02
termination code = 0
DIMACS: 1.4e-13 0.0e+00 8.0e-08 0.0e+00 3.2e-09 3.8e-09
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +59.6955

Figure 1 Histogram of residual amplitudes for ℓ1-norm, with the scaled penalty functions

2) For ℓ2-norm, the object function is Ax-b2 , and using the following code to solve the problem
cvx_begin
variable x(n);
minimize( norm(A*x-b,2));
cvx_end
The running result is as following:
Calling SDPT3 4.0: 101 variables, 31 equality constraints
For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

num. of constraints = 31
dim. of socp var = 101, num. of socp blk = 1
*******************************************************************
SDPT3: Infeasible path-following algorithms
*******************************************************************
version predcorr gam expon scale_data
NT 1 0.000 1 0
it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime
-------------------------------------------------------------------
0|0.000|0.000|4.5e+00|1.4e+00|1.2e+02| 0.000000e+00 0.000000e+00| 0:0:00| chol 1 1
1|0.916|1.000|3.8e-01|8.4e-03|1.6e+01|-1.107217e+01 -1.558450e+01| 0:0:00| chol 1 1
2|1.000|1.000|7.7e-08|8.4e-04|1.3e+00|-8.059478e+00 -9.392154e+00| 0:0:00| chol 1 1
3|0.989|0.986|1.9e-08|9.5e-05|1.7e-02|-8.467255e+00 -8.483532e+00| 0:0:00| chol 1 1
4|0.989|0.989|5.8e-09|9.4e-06|1.9e-04|-8.472550e+00 -8.472631e+00| 0:0:00| chol 1 1
5|0.989|0.989|6.4e-11|1.0e-07|2.1e-06|-8.472608e+00 -8.472609e+00| 0:0:00| chol 1 1
6|0.989|0.989|7.0e-13|1.1e-09|2.4e-08|-8.472609e+00 -8.472609e+00| 0:0:00|
stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
number of iterations = 6
primal objective value = -8.47260907e+00
dual objective value = -8.47260908e+00
gap := trace(XZ) = 2.35e-08
relative gap = 1.31e-09
actual relative gap = 5.50e-10
rel. primal infeas = 7.01e-13
rel. dual infeas = 1.15e-09
norm(X), norm(y), norm(Z) = 1.4e+00, 8.5e+00, 1.2e+01
norm(A), norm(b), norm(C) = 5.5e+01, 2.0e+00, 1.2e+01
Total CPU time (secs) = 0.08
CPU time per iteration = 0.01
termination code = 0
DIMACS: 7.0e-13 0.0e+00 3.9e-09 0.0e+00 5.5e-10 1.3e-09
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +8.47261

Figure 2 Histogram of residual amplitudes for ℓ2-norm, with the (scaled) penalty function
3) For deadzone linear, the object function is sum(ϕu),where
ϕu=0, &u≤0.5u-0.5, &u>0.5
u=Ax-b and using the following code to solve the problem
cvx_begin
variable x3(n);
minimize( sum(max(abs(A*x3-b)-0.5,0)));
cvx_end
The running result is as following:
Calling SDPT3 4.0: 430 variables, 200 equality constraints
------------------------------------------------------------

num. of constraints = 200
dim. of socp var = 200, num. of socp blk = 100
dim. of linear var = 200
dim. of free var = 30 *** convert ublk to lblk
*******************************************************************
SDPT3: Infeasible path-following algorithms
*******************************************************************
version predcorr gam expon scale_data
NT 1 0.000 1 0
it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime
-------------------------------------------------------------------
0|0.000|0.000|1.8e+00|6.1e+01|5.9e+05| 4.994495e+03 0.000000e+00| 0:0:00| chol 1 1
1|1.000|0.430|4.1e-05|3.5e+01|1.8e+05| 1.298792e+04 -3.753191e+02| 0:0:00| chol 1 1
2|1.000|0.986|2.9e-05|5.6e-01|1.3e+04| 1.027122e+04 -5.033511e+01| 0:0:00| chol 1 1
3|0.998|1.000|1.5e-06|2.7e-02|3.2e+02| 2.672257e+02 -4.306854e+01| 0:0:00| chol 1 1
4|0.929|0.692|1.2e-05|1.0e-02|8.7e+01| 6.834936e+01 -1.649617e+01| 0:0:00| chol 1 1
5|0.932|0.267|4.9e-06|7.5e-03|5.2e+01| 4.315547e+01 -7.356107e+00| 0:0:00| chol 1 1
6|0.971|0.389|3.1e-06|4.6e-03|3.7e+01| 4.304656e+01 7.233939e+00| 0:0:00| chol 1 1
7|1.000|0.319|4.9e-07|3.1e-03|2.1e+01| 3.319546e+01 1.265459e+01| 0:0:00| chol 1 1
8|1.000|0.482|1.3e-07|1.6e-03|1.1e+01| 2.949734e+01 1.903483e+01| 0:0:00| chol 1 1
9|1.000|0.477|3.9e-08|8.5e-04|5.2e+00| 2.779047e+01 2.265770e+01| 0:0:00| chol 1 1
10|1.000|0.287|2.8e-08|6.1e-04|3.9e+00| 2.758221e+01 2.378289e+01| 0:0:00| chol 1 1
11|0.978|0.427|1.5e-08|3.5e-04|2.2e+00| 2.726306e+01 2.507759e+01| 0:0:00| chol 1 1
12|1.000|0.467|6.4e-09|1.9e-04|1.2e+00| 2.710945e+01 2.594507e+01| 0:0:00| chol 1 1
13|1.000|0.497|1.6e-09|9.3e-05|5.8e-01| 2.704017e+01 2.646287e+01| 0:0:00| chol 1 1
14|1.000|0.167|2.5e-10|7.8e-05|5.2e-01| 2.706103e+01 2.654623e+01| 0:0:00| chol 1 1
15|0.946|0.485|5.9e-10|4.0e-05|2.7e-01| 2.702852e+01 2.676251e+01| 0:0:00| chol 1 1
16|1.000|0.342|1.5e-09|2.6e-05|1.8e-01| 2.702033e+01 2.684315e+01| 0:0:00| chol 1 1
17|1.000|0.554|4.1e-10|1.2e-05|7.9e-02| 2.701150e+01 2.693328e+01| 0:0:00| chol 1 1
18|1.000|0.376|2.6e-11|7.4e-06|5.0e-02| 2.701008e+01 2.696100e+01| 0:0:00| chol 1 1
19|1.000|0.430|2.7e-10|4.2e-06|2.8e-02| 2.700907e+01 2.698114e+01| 0:0:00| chol 1 1
20|0.972|0.824|2.7e-11|2.2e-05|5.2e-03| 2.700845e+01 2.700357e+01| 0:0:00| chol 1 1
21|0.987|0.985|3.6e-13|4.0e-06|9.6e-05| 2.700839e+01 2.700832e+01| 0:0:00| chol 1 1
22|1.000|0.989|3.7e-15|7.3e-08|1.9e-06| 2.700839e+01 2.700839e+01| 0:0:00| chol 1 1
23|1.000|0.989|6.9e-15|1.4e-09|3.2e-08| 2.700839e+01 2.700839e+01| 0:0:00|
stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
number of iterations = 23
primal objective value = 2.70083890e+01
dual objective value = 2.70083889e+01
gap := trace(XZ) = 3.24e-08
relative gap = 5.89e-10
actual relative gap = 4.61e-10
rel. primal infeas = 6.88e-15
rel. dual infeas = 1.42e-09
norm(X), norm(y), norm(Z) = 1.4e+01, 9.3e+00, 1.3e+01
norm(A), norm(b), norm(C) = 8.0e+01, 1.3e+01, 1.1e+01
Total CPU time (secs) = 0.43
CPU time per iteration = 0.02
termination code = 0
DIMACS: 2.5e-14 0.0e+00 7.8e-09 0.0e+00 4.6e-10 5.9e-10
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +27.0084

Figure 3 Histogram of residual amplitudes for deadzone-linear, with the (scaled) penalty function
4) For log-barrier, the object function is sum(ϕu), where
ϕu=-log⁡(1-u2), &u<1∞ , &u≥1
and u=Ax-b and using the following code to solve the problem

cvx_begin
variable x4(n);
u = A*x-b;
minimize( sum(-1*log(1-u.^2)));
subject to
abs(u)<1
cvx_end

The running result is as following:
Successive approximation method to be employed.
For improved efficiency, SDPT3 is solving the dual problem.
SDPT3 will be called several times to refine the solution.
Original size: 900 variables, 330 equality constraints
100 exponentials add 800 variables, 500 equality constraints
-----------------------------------------------------------------
Cones | Errors |
Mov/Act | Centering Exp cone Poly cone | Status
--------+---------------------------------+---------
100/100 | 3.250e+00 6.195e-01 0.000e+00 | Solved
100/100 | 6.211e-01 2.532e-02 0.000e+00 | Solved
100/100 | 5.379e-02 1.879e-04 0.000e+00 | Solved
100/100 | 6.468e-03 2.721e-06 0.000e+00 | Solved
24/100 | 7.769e-04 3.731e-08 0.000e+00 | Solved
0/ 0 | 0.000e+00 0.000e+00 0.000e+00 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +33.9336

Figure 4 Histogram of residual amplitudes for log-barrier, with the (scaled) penalty function, and the quadratic penalty is also shown, in dashed curve.
Conclusion
From the plots of the penalty functions we note that
l The ℓ1-norm penalty puts the most weight on small residuals and the least weight on large residuals.
l The ℓ2-norm penalty puts very small weight on small residuals, but strong weight on large residuals.
l The deadzone-linear penalty function puts no weight on residuals smaller than 0.5, and relatively little weight on large residuals.
l The log barrier penalty puts weight very much like the ℓ2-norm penalty for small residuals, but puts very strong weight on residuals larger than around 0.8, and infinite weight on residuals larger than 1.
Several features are clear from the amplitude distributions:
l For the ℓ1-optimal solution, many residuals are either zero or very small. The ℓ1-optimal solution also has relatively more large residuals.
l The ℓ2-norm approximation has many modest residuals, and relatively few larger ones.
l For the deadzone-linear penalty, we see that many residuals have the value ±0.5, right at the edge of the ‘free’ zone, for which no penalty is assessed.
l For the log barrier penalty, we see that no residuals have a magnitude larger than 1, but otherwise the residual distribution is similar to the residual distribution for ℓ2-norm approximation.
According the content at Page 298, I plot figure 6.3 as following:

Figure 5 A (nonconvex) penalty function
Example 6.3 Optimal input design
We consider a dynamical system with scalar input sequence u0,u1,…,u(N), and scalar output sequence y0,y1,…,y(N), related by convolution:
yt=τ=0tuτht-τ,t=0,1,2,…,N.
The sequence h(0), h(1),… , h(N) is called the convolution kernel or impulse response of the system.
Our goal is to choose the input sequence u to achieve several goals.
l Output tracking. The primary goal is that the output y should track, or follow, a desired target or reference signal ydes. We measure output tracking error by the quadratic function
Jtrack=1N+1t=0N[yt-ydes(t)]2
l Smallest input. The input should not be large. We measure the magnitude of the input by the quadratic function
Jmag=1N+1t=0N[ut]2
l Smallest input variations. The input should not vary rapidly. We measure the magnitude of the input variations by the quadratic function
Jder=1N+1t=0N-1[ut+1-u(t)]2
By minimizing a weighted sum
Jtrack+δJder+ηJmag
where δ>0 and η>0, we can trade off the three objectives.
Now