9.4.2 多目的制御


 アームシステムの多目的制御 (ソルバ SeDuMi を利用) ...... M ファイル LMI_multi_objective.m (p.195)
clear
format compact; format short e
% ----------
J  = 0.0712;  M  = 0.390;  mu = 0.695;
l  = 0.204;   g  = 9.81;

A = [     0      1
      M*g*l/J  -mu/J ];
B = [  0
      1/J ];
% ----------
n = 2;  p = 1;
X = sdpvar(n,n,'sy');
F = sdpvar(p,n,'f');
Z = sdpvar(n,n,'sy');
% ----------
M = A*X + B*F;
LMI = set([]);
% ----------
qc = -5;  rc =  4;

LMI = LMI + set([qc*(M+M')+(rc^2-qc^2)*X  M
                           M'             X ] > 0);
% ----------
Qh = diag([sqrt(10) 0]);  R  = 1;

LMI = LMI + set([-(M+M')    X*Qh        F'*R
                  Qh*X     eye(n)    zeros(n,p)
                   R*F   zeros(p,n)      R     ] > 0);
LMI = LMI + set([   Z     eye(n)
                  eye(n)    X   ] > 0);
% ----------
options = sdpsettings('solver','sedumi');
solvesdp(LMI,trace(Z),options) 

Z_opt = double(Z); trace(Z_opt)
X_opt = double(X)
F_opt = double(F)
K_opt = F_opt*inv(X_opt)
% ----------
format short
>> LMI_multi_objective
SeDuMi 1.3 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, theta = 0.250, beta = 0.500
eqs m = 8, order n = 14, dim = 58, blocks = 4
nnz(A) = 38 + 0, nnz(ADA) = 52, nnz(L) = 30
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            1.47E+003 0.000
  1 : -1.16E+000 4.46E+002 0.000 0.3042 0.9000 0.9000   1.63  1  1  8.7E+001
  2 : -3.04E+000 1.34E+002 0.000 0.2996 0.9000 0.9000   0.50  1  1  3.9E+001
  3 : -4.61E+000 3.75E+001 0.000 0.2807 0.9000 0.9000   0.49  1  1  1.5E+001
  4 : -5.61E+000 1.03E+001 0.000 0.2738 0.9000 0.9000   0.41  1  1  6.1E+000
  5 : -5.95E+000 3.10E+000 0.000 0.3019 0.9000 0.9000   0.44  1  1  2.5E+000
  6 : -5.82E+000 9.43E-001 0.000 0.3039 0.9000 0.9000   0.46  1  1  1.1E+000
  7 : -5.53E+000 3.15E-001 0.000 0.3337 0.9000 0.9000   0.37  1  1  5.3E-001
  8 : -5.36E+000 9.93E-002 0.000 0.3156 0.9000 0.9000   0.45  1  1  2.3E-001
  9 : -5.20E+000 3.70E-002 0.000 0.3728 0.9000 0.9000   0.43  1  1  1.2E-001
 10 : -5.11E+000 1.27E-002 0.000 0.3436 0.9000 0.9000   0.55  1  1  5.3E-002
 11 : -5.05E+000 4.32E-003 0.000 0.3392 0.9000 0.9000   0.60  1  1  2.2E-002
 12 : -5.03E+000 9.63E-004 0.000 0.2231 0.9000 0.9000   0.76  1  1  5.7E-003
 13 : -5.02E+000 5.65E-005 0.000 0.0587 0.9900 0.9900   0.90  1  1  3.5E-004
 14 : -5.02E+000 4.14E-006 0.404 0.0733 0.9900 0.9900   0.99  1  1  2.6E-005
 15 : -5.02E+000 1.01E-006 0.000 0.2436 0.9000 0.9000   1.00  1  1  6.4E-006
 16 : -5.02E+000 5.40E-008 0.000 0.0535 0.9900 0.9900   1.00  1  1  3.4E-007
 17 : -5.02E+000 1.73E-008 0.000 0.3204 0.9000 0.9000   1.00  2  2  1.1E-007
 18 : -5.02E+000 3.33E-009 0.000 0.1925 0.9000 0.9000   1.00  2  2  2.1E-008
 19 : -5.02E+000 6.57E-010 0.000 0.1974 0.9000 0.9000   1.00  2  2  4.1E-009
 20 : -5.02E+000 1.55E-010 0.000 0.2364 0.9000 0.9000   1.00  2  2  9.8E-010

iter seconds digits       c*x               b*y
 20      0.4  10.0 -5.0185569346e+000 -5.0185569341e+000
|Ax-b| =  1.4e-009, [Ay-c]_+ =  0.0E+000, |x|= 2.7e+001, |y|= 1.7e+001

Detailed timing (sec)
   Pre          IPM          Post
7.001E-003    2.650E-001    2.997E-003    
Max-norms: ||b||=1, ||c|| = 2,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 185619.
ans = 
    yalmiptime: 1.1700e-001
    solvertime: 2.7600e-001
          info: 'No problems detected (SeDuMi-1.3)'
       problem: 0
        dimacs: [7.1804e-010 0 0 0 -4.3415e-011 5.2496e-010]
ans =
  5.0186e+000
X_opt =
  4.5315e-001 -1.9664e+000
 -1.9664e+000  1.5586e+001
F_opt =
 -1.1853e+000  5.4031e+000
K_opt =
 -2.4558e+000  3.6822e-002

 アームシステムの多目的制御 (ソルバ SDPT3 を利用)
   ...... M ファイル LMI_multi_objective_sdpt3.m(本には記載していない)
clear
format compact; format short e
% ----------
J  = 0.0712;  M  = 0.390;  mu = 0.695;
l  = 0.204;   g  = 9.81;

A = [     0      1
      M*g*l/J  -mu/J ];
B = [  0
      1/J ];
% ----------
n = 2;  p = 1;
X = sdpvar(n,n,'sy');
F = sdpvar(p,n,'f');
Z = sdpvar(n,n,'sy');
% ----------
M = A*X + B*F;
LMI = set([]);
% ----------
qc = -5;  rc =  4;

LMI = LMI + set([qc*(M+M')+(rc^2-qc^2)*X  M
                           M'             X ] > 0);
% ----------
Qh = diag([sqrt(10) 0]);  R  = 1;

LMI = LMI + set([-(M+M')    X*Qh        F'*R
                  Qh*X     eye(n)    zeros(n,p)
                   R*F   zeros(p,n)      R     ] > 0);
LMI = LMI + set([   Z     eye(n)
                  eye(n)    X   ] > 0);
% ----------
options = sdpsettings('solver','sdpt3');
solvesdp(LMI,trace(Z),options) 

Z_opt = double(Z); trace(Z_opt)
X_opt = double(X)
F_opt = double(F)
K_opt = F_opt*inv(X_opt)
% ----------
format short
>> LMI_multi_objective_sdpt3

 num. of constraints =  8
 dim. of sdp    var  = 13,   num. of sdp  blk  =  3
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
   HKM      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
 0|0.000|0.000|1.0e+003|8.0e+001|7.6e+003| 3.000000e+001  0.000000e+000| 0:0:00| chol  1  1 
 1|0.749|0.723|2.6e+002|2.2e+001|2.5e+003| 4.237109e+001 -3.653107e+001| 0:0:00| chol  1  1 
 2|0.744|0.724|6.6e+001|6.2e+000|9.7e+002| 7.645543e+001 -6.305254e+001| 0:0:00| chol  1  1 
 3|0.808|0.799|1.3e+001|1.2e+000|2.6e+002| 5.464762e+001 -3.830663e+001| 0:0:00| chol  1  1 
 4|0.712|0.935|3.6e+000|8.1e-002|8.9e+001| 3.787390e+001 -2.134200e+001| 0:0:00| chol  1  1 
 5|0.848|1.000|5.6e-001|9.9e-006|2.8e+001| 1.198935e+001 -1.086914e+001| 0:0:00| chol  1  1 
 6|0.827|1.000|9.6e-002|9.9e-007|6.2e+000|-9.884840e-001 -6.223500e+000| 0:0:00| chol  1  1 
 7|0.897|1.000|9.9e-003|4.5e-003|1.9e+000|-3.913768e+000 -5.507548e+000| 0:0:00| chol  1  1 
 8|0.949|0.987|5.0e-004|3.8e-004|1.2e-001|-4.934380e+000 -5.039692e+000| 0:0:00| chol  1  1 
 9|0.976|0.974|1.2e-005|2.0e-005|3.9e-003|-5.016152e+000 -5.019185e+000| 0:0:00| chol  1  1 
10|0.977|0.980|2.8e-007|6.3e-007|9.3e-005|-5.018504e+000 -5.018572e+000| 0:0:00| chol  1  1 
11|0.921|0.974|2.2e-008|3.3e-008|6.4e-006|-5.018553e+000 -5.018558e+000| 0:0:00| chol  1  2 
12|1.000|1.000|1.1e-012|3.4e-009|1.3e-006|-5.018556e+000 -5.018557e+000| 0:0:00| chol  1  1 
13|1.000|1.000|1.1e-010|1.0e-012|7.0e-008|-5.018557e+000 -5.018557e+000| 0:0:00|
  stop: max(relative gap, infeasibilities) < 1.00e-007
-------------------------------------------------------------------
 number of iterations   = 13
 primal objective value = -5.01855689e+000
 dual   objective value = -5.01855696e+000
 gap := trace(XZ)       = 7.04e-008
 relative gap           = 6.38e-009
 actual relative gap    = 6.28e-009
 rel. primal infeas     = 1.08e-010
 rel. dual   infeas     = 1.00e-012
 norm(X), norm(y), norm(Z) = 2.7e+001, 1.7e+001, 8.8e+002
 norm(A), norm(b), norm(C) = 2.5e+002, 2.4e+000, 3.6e+000
 Total CPU time (secs)  = 0.4  
 CPU time per iteration = 0.0  
 termination code       =  0
 DIMACS: 1.3e-010  0.0e+000  1.8e-012  0.0e+000  6.3e-009  6.4e-009
-------------------------------------------------------------------
ans = 
    yalmiptime: 1.1500e-001
    solvertime: 5.0200e-001
          info: 'No problems detected (SDPT3-4)'
       problem: 0
        dimacs: [1.3036e-010 0 3.1583e-012 0 6.2790e-009 6.3728e-009]
ans =
  5.0186e+000
X_opt =
  4.5315e-001 -1.9665e+000
 -1.9665e+000  1.5586e+001
F_opt =
 -1.1853e+000  5.4031e+000
K_opt =
 -2.4559e+000  3.6812e-002

 アームシステムの多目的制御 (Robust Control Toolbox のソルバ LMILAB を利用)
   ...... M ファイル LMI_multi_objective_lmilab.m(本には記載していない)
clear
format compact; format short e
% ----------
J  = 0.0712;  M  = 0.390;  mu = 0.695;
l  = 0.204;   g  = 9.81;

A = [     0      1
      M*g*l/J  -mu/J ];
B = [  0
      1/J ];
% ----------
n = 2;  p = 1;
X = sdpvar(n,n,'sy');
F = sdpvar(p,n,'f');
Z = sdpvar(n,n,'sy');
% ----------
M = A*X + B*F;
LMI = set([]);
% ----------
qc = -5;  rc =  4;

LMI = LMI + set([qc*(M+M')+(rc^2-qc^2)*X  M
                           M'             X ] > 0);
% ----------
Qh = diag([sqrt(10) 0]);  R  = 1;

LMI = LMI + set([-(M+M')    X*Qh        F'*R
                  Qh*X     eye(n)    zeros(n,p)
                   R*F   zeros(p,n)      R     ] > 0);
LMI = LMI + set([   Z     eye(n)
                  eye(n)    X   ] > 0);
% ----------
options = sdpsettings('solver','lmilab');
solvesdp(LMI,trace(Z),options) 

Z_opt = double(Z); trace(Z_opt)
X_opt = double(X)
F_opt = double(F)
K_opt = F_opt*inv(X_opt)
% ----------
format short
>> LMI_multi_objective_lmilab

 Solver for linear objective minimization under LMI constraints 

 Iterations   :    Best objective value so far 
 
     1
     2
     3
     4
     5                  15.917321
     6                  11.322246
     7                   8.938944
     8                   8.316948
     9                   7.321659
    10                   7.321659
    11                   7.321659
    12                   5.705634
    13                   5.705634
    14                   5.234645
***                 new lower bound:     3.673585
    15                   5.205537
***                 new lower bound:     4.224103
    16                   5.163084
***                 new lower bound:     4.562963
    17                   5.131253
***                 new lower bound:     4.729848
    18                   5.096550
***                 new lower bound:     4.834394
    19                   5.096550
***                 new lower bound:     4.900231
    20                   5.036927
***                 new lower bound:     4.994248
    21                   5.024449
***                 new lower bound:     4.998723
    22                   5.022059
***                 new lower bound:     5.006302
    23                   5.020604
***                 new lower bound:     5.010957
    24                   5.019930
***                 new lower bound:     5.015161

 Result:  feasible solution of required accuracy
          best objective value:     5.019930
          guaranteed relative accuracy: 9.50e-004
          f-radius saturation:  0.000% of R = 1.00e+009 
 
ans = 
    yalmiptime: 9.1800e-001
    solvertime: 4.2000e-001
          info: 'No problems detected (LMILAB)'
       problem: 0
        dimacs: [NaN NaN 0 0 NaN NaN]
ans =
  5.0199e+000
X_opt =
  4.5313e-001 -1.9663e+000
 -1.9663e+000  1.5583e+001
F_opt =
 -1.1856e+000  5.4043e+000
K_opt =
 -2.4565e+000  3.6838e-002

 例 9.3 (p.188):LMI による円領域への極配置 (ソルバ SeDuMi を利用)
   ...... M ファイル LMI_poleplacement.m (p.196 の下 3 行目)
clear
format compact; format short e
% ----------
J  = 0.0712;  M  = 0.390;  mu = 0.695;
l  = 0.204;   g  = 9.81;

A = [     0      1
      M*g*l/J  -mu/J ];
B = [  0
      1/J ];
% ----------
n = 2;  p = 1;
X = sdpvar(n,n,'sy');
F = sdpvar(p,n,'f');
% ----------
M = A*X + B*F;
LMI = set([]);
% ----------
qc = -5;  rc =  4;

LMI = LMI + set([qc*(M+M')+(rc^2-qc^2)*X  M
                           M'             X ] > 0);
% ----------
options = sdpsettings('solver','sedumi');
solvesdp(LMI,[],options) 

X_opt = double(X)
F_opt = double(F)
K_opt = F_opt*inv(X_opt)
% ----------
format short
>> LMI_poleplacement
SeDuMi 1.3 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, theta = 0.250, beta = 0.500
eqs m = 5, order n = 5, dim = 17, blocks = 2
nnz(A) = 20 + 0, nnz(ADA) = 25, nnz(L) = 15
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            1.61E+004 0.000
  1 :  0.00E+000 2.94E+003 0.000 0.1821 0.9000 0.9000   1.00  1  0  5.2E+001
  2 :  0.00E+000 2.72E+002 0.000 0.0927 0.9900 0.9900   1.00  1  1  4.8E+000
  3 :  0.00E+000 1.97E+001 0.000 0.0723 0.9900 0.9900   1.00  1  1  3.5E-001
  4 :  0.00E+000 5.36E-003 0.000 0.0003 0.9999 0.9999   1.00  1  1  9.4E-005
  5 :  0.00E+000 7.80E-010 0.000 0.0000 1.0000 1.0000   1.00  1  1  2.9E-011

iter seconds digits       c*x               b*y
  5      0.1   Inf  0.0000000000e+000  0.0000000000e+000
|Ax-b| =  9.9e-012, [Ay-c]_+ =  0.0E+000, |x|= 2.9e-012, |y|= 1.1e+000

Detailed timing (sec)
   Pre          IPM          Post
7.001E-003    5.300E-002    2.002E-003    
Max-norms: ||b||=0, ||c|| = 0,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.47928.
ans = 
    yalmiptime: 1.1600e-001
    solvertime: 6.3000e-002
          info: 'No problems detected (SeDuMi-1.3)'
       problem: 0
        dimacs: [9.8964e-012 0 0 0 0 2.0854e-013]
X_opt =
  3.2038e-001 -3.9160e-001
 -3.9160e-001  5.9828e-001
F_opt =
 -5.0256e-001  6.6852e-001
K_opt =
 -1.0144e+000  4.5342e-001

 例 9.4 (p.191):LMI による最適レギュレータ設計 (ソルバ SeDuMi を利用)
   ...... M ファイル LMI_LQ.m (p.196 の下 1 行目)
clear
format compact; format short e
% ----------
J  = 0.0712;  M  = 0.390;  mu = 0.695;
l  = 0.204;   g  = 9.81;

A = [     0      1
      M*g*l/J  -mu/J ];
B = [  0
      1/J ];
% ----------
n = 2;  p = 1;
X = sdpvar(n,n,'sy');
F = sdpvar(p,n,'f');
% ----------
M = A*X + B*F;
LMI = set([]);
% ----------
qc = -5;  rc =  4;

LMI = LMI + set([qc*(M+M')+(rc^2-qc^2)*X  M
                           M'             X ] > 0);
% ----------
options = sdpsettings('solver','sedumi');
solvesdp(LMI,[],options) 

X_opt = double(X)
F_opt = double(F)
K_opt = F_opt*inv(X_opt)
% ----------
format short
>> LMI_LQ
SeDuMi 1.3 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, theta = 0.250, beta = 0.500
eqs m = 8, order n = 10, dim = 42, blocks = 3
nnz(A) = 18 + 0, nnz(ADA) = 52, nnz(L) = 30
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            8.21E+001 0.000
  1 : -1.09E+000 2.49E+001 0.000 0.3032 0.9000 0.9000   1.63  1  1  1.5E+001
  2 : -2.78E+000 6.70E+000 0.000 0.2693 0.9000 0.9000   0.55  1  1  5.8E+000
  3 : -4.03E+000 1.73E+000 0.000 0.2577 0.9000 0.9000   0.51  1  1  2.0E+000
  4 : -4.66E+000 4.57E-001 0.000 0.2646 0.9000 0.9000   0.46  1  1  7.7E-001
  5 : -4.72E+000 1.34E-001 0.000 0.2924 0.9000 0.9000   0.45  1  1  3.2E-001
  6 : -4.58E+000 4.47E-002 0.000 0.3349 0.9000 0.9000   0.46  1  1  1.5E-001
  7 : -4.47E+000 1.55E-002 0.000 0.3470 0.9000 0.9000   0.51  1  1  6.8E-002
  8 : -4.33E+000 5.91E-003 0.000 0.3809 0.9000 0.9000   0.41  1  1  3.9E-002
  9 : -4.26E+000 2.16E-003 0.000 0.3650 0.9000 0.9000   0.53  1  1  1.8E-002
 10 : -4.16E+000 7.78E-004 0.000 0.3604 0.9000 0.9000   0.39  1  1  1.0E-002
 11 : -4.11E+000 2.78E-004 0.000 0.3579 0.9000 0.9000   0.49  1  1  4.8E-003
 12 : -4.05E+000 9.77E-005 0.000 0.3509 0.9000 0.9000   0.35  1  1  2.7E-003
 13 : -4.01E+000 3.36E-005 0.000 0.3435 0.9000 0.9000   0.47  1  1  1.2E-003
 14 : -3.98E+000 1.10E-005 0.000 0.3269 0.9000 0.9000   0.35  1  1  6.5E-004
 15 : -3.96E+000 3.74E-006 0.000 0.3411 0.9000 0.9000   0.46  1  1  3.0E-004
 16 : -3.94E+000 1.18E-006 0.000 0.3152 0.9000 0.9000   0.35  1  1  1.5E-004
 17 : -3.93E+000 3.68E-007 0.000 0.3118 0.9000 0.9000   0.51  1  1  6.2E-005
 18 : -3.92E+000 1.42E-007 0.000 0.3856 0.9000 0.9000   0.53  1  1  3.0E-005
 19 : -3.92E+000 4.91E-009 0.000 0.0346 0.9900 0.9900   0.90  1  1  1.1E-006
 20 : -3.92E+000 2.43E-010 0.308 0.0496 0.9900 0.9900   0.99  2  2  5.3E-008
 21 : -3.92E+000 9.02E-012 0.000 0.0371 0.9900 0.9900   1.00  2  2  2.0E-009
 22 : -3.92E+000 7.89E-013 0.152 0.0875 0.9900 0.9900   1.00  2  2  1.7E-010

iter seconds digits       c*x               b*y
 22      0.5   9.1 -3.9164628691e+000 -3.9164628659e+000
|Ax-b| =  2.5e-010, [Ay-c]_+ =  0.0E+000, |x|= 1.7e+001, |y|= 4.0e+002

Detailed timing (sec)
   Pre          IPM          Post
7.001E-003    2.600E-001    2.002E-003    
Max-norms: ||b||=1, ||c|| = 2,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 275131.
ans = 
    yalmiptime: 1.1600e-001
    solvertime: 2.7000e-001
          info: 'No problems detected (SeDuMi-1.3)'
       problem: 0
        dimacs: [1.2571e-010 0 0 0 -3.6945e-010 5.8581e-009]
ans =
  3.9165e+000
X_opt =
  2.4208e+000 -2.9301e+001
 -2.9301e+001  3.9675e+002
F_opt =
 -2.8247e-004 -1.4041e+001
K_opt =
 -4.0377e+000 -3.3358e-001