Merge pull request #895 from jrforbes/main · python-control/python-control@fe4589b

1+

"""H2 design using h2syn.

2+3+

Demonstrate H2 desing for a SISO plant using h2syn. Based on [1], Ex. 7.

4+5+

[1] Scherer, Gahinet, & Chilali, "Multiobjective Output-Feedback Control via

6+

LMI Optimization", IEEE Trans. Automatic Control, Vol. 42, No. 7, July 1997.

7+8+

[2] Zhou & Doyle, "Essentials of Robust Control", Prentice Hall,

9+

Upper Saddle River, NJ, 1998.

10+

"""

11+

# %%

12+

# Packages

13+

import numpy as np

14+

import control

15+16+

# %%

17+

# State-space system.

18+19+

# Process model.

20+

A = np.array([[0, 10, 2],

21+

[-1, 1, 0],

22+

[0, 2, -5]])

23+

B1 = np.array([[1],

24+

[0],

25+

[1]])

26+

B2 = np.array([[0],

27+

[1],

28+

[0]])

29+30+

# Plant output.

31+

C2 = np.array([[0, 1, 0]])

32+

D21 = np.array([[2]])

33+

D22 = np.array([[0]])

34+35+

# H2 performance.

36+

C1 = np.array([[0, 1, 0],

37+

[0, 0, 1],

38+

[0, 0, 0]])

39+

D11 = np.array([[0],

40+

[0],

41+

[0]])

42+

D12 = np.array([[0],

43+

[0],

44+

[1]])

45+46+

# Dimensions.

47+

n_u, n_y = 1, 1

48+49+

# %%

50+

# H2 design using h2syn.

51+52+

# Create augmented plant.

53+

Baug = np.block([B1, B2])

54+

Caug = np.block([[C1], [C2]])

55+

Daug = np.block([[D11, D12], [D21, D22]])

56+

Paug = control.ss(A, Baug, Caug, Daug)

57+58+

# Input to h2syn is Paug, number of inputs to controller,

59+

# and number of outputs from the controller.

60+

K = control.h2syn(Paug, n_y, n_u)

61+62+

# Extarct controller ss realization.

63+

A_K, B_K, C_K, D_K = K.A, K.B, K.C, K.D

64+65+

# %%

66+

# Compute closed-loop H2 norm.

67+68+

# Compute closed-loop system, Tzw(s). See Eq. 4 in [1].

69+

Azw = np.block([[A + B2 @ D_K @ C2, B2 @ C_K],

70+

[B_K @ C2, A_K]])

71+

Bzw = np.block([[B1 + B2 @ D_K @ D21],

72+

[B_K @ D21]])

73+

Czw = np.block([C1 + D12 @ D_K @ C2, D12 @ C_K])

74+

Dzw = D11 + D12 @ D_K @ D21

75+

Tzw = control.ss(Azw, Bzw, Czw, Dzw)

76+77+

# Compute closed-loop H2 norm via Lyapunov equation.

78+

# See [2], Lemma 4.4, pg 53.

79+

Qzw = control.lyap(Azw.T, Czw.T @ Czw)

80+

nu = np.sqrt(np.trace(Bzw.T @ Qzw @ Bzw))

81+

print(f'The closed-loop H_2 norm of Tzw(s) is {nu}.')

82+

# Value is 7.748350599360575, the same as reported in [1].

83+84+

# %%