# How do I design a close-loop controller on my simulink of piezo bender of energy harvester

3 次查看（过去 30 天）
Fajwa Noor 2022-6-14

The aim is to design a more sophisticated close-loop controller to optimize the transfer of power and improve the efficiency of the energy harvester in different conditions.
##### 4 个评论显示 2更早的评论隐藏 2更早的评论
Sam Chak 2022-6-16
Thanks. Found it. Now, what type of "sophistication" that you want to design for ?
Fajwa Noor 2022-6-16
I want to design a close-loop controller that can optimize the transfer of power and improve the efficiency of the energy harvester if the vibration source does not have a constant frequency

### 采纳的回答

Sam Chak 2022-6-16

Generally, the optimization should come later once you have successfully designed the basic structure of the mathematical law that stabilizes the closed-loop system. If the state-feedback mathematical law is acceptable to you, then you may consider designing multiple Linear Quadratic Regulators (LQRs) for the piezobender when the vibration source have different frequencies.
Basic info about using the K = LQR(A, B, Q, R) command in MATLAB can be found here:
See the example of designing a family of PID controllers for multiple operating points (which I think it is similar to your case of having a vibration source with different frequencies)
EDIT: To include some educational materials.
Part 1: Mathematical preliminaries
From the dynamical equations of the piezobender, we want to design {} to stabilize the system
But first, we need to rewrite the equations in a compact form called the matrix differential equation:
where the vector and the generalized force vector .
Rearranging the matrix differential equation yields
.
Now, we want to express the matrix differential equation in state-space representation so that the LQR can be designed in MATLAB as described above
.
Note that is chosen as the math notation for the identity matrix.
The problem is reformulated to designing a state-feedback mathematical law given by
,
such that
is a Hurwitz matrix, having the characteristic polynomial with all its roots have strictly negative real part. The design process can be fairly tedious in pen and paper, noting that the size of the state matrix is . So, in such situation, the LQR algorithm can be used to calculate the optimal Gain matrix so that the following quadratic cost function
is minimized. The LQR algorithm requires and to be defined. As a starter (to see how well the LQR performs), we can select
... matching the size of matrix
... due to has 4 inputs.
Part 2: LQR design
% Parameters
rho = 7804.89;
l = 3.175e-2;
w = 1.72e-2;
d = 5.084e-4;
E = 1.36628e11;
I = 1.3907e-14;
e31 = 7.5459;
epsl = 1.38965e-8;
bm = 0;
kb = 1e-5;
% Mass matrix
m = (rho*l*w*d)/420;
M1 = [156 22*l 54 -13*l;
22*l 4*l^2 13*l -3*l^2;
54 13*l 156 -22*l;
-13*l -3*l^2 -22*l 4*l^2];
M = m*M1;
% Stiffness matrix
K = [12*E*I/l^3 6*E*I/l^2 -12*E*I/l^3 6*E*I/l^2;
6*E*I/l^2 4*E*I/l -6*E*I/l^2 2*E*I/l;
-12*E*I/l^3 -6*E*I/l^2 12*E*I/l^3 -6*E*I/l^2;
6*E*I/l^2 2*E*I/l -6*E*I/l^2 4*E*I/l];
% Damping matrix
B = kb*K + bm*M;
State-space
n = size(M, 1);
a11 = zeros(n);
a12 = eye(n);
a21 = - M\K;
a22 = - M\B;
b1 = zeros(n);
b2 = M\a12;
A = [a11 a12; a21 a22];
B = [b1; b2];
C = [eye(n) zeros(n)];
D = zeros(n);
sys1 = ss(A, B, C, D)
sys1 = A = x1 x2 x3 x4 x5 x6 x7 x8 x1 0 0 0 0 1 0 0 0 x2 0 0 0 0 0 1 0 0 x3 0 0 0 0 0 0 1 0 x4 0 0 0 0 0 0 0 1 x5 2.301e+07 4.175e+05 -2.301e+07 3.131e+05 230.1 4.175 -230.1 3.131 x6 -8.698e+09 -1.479e+08 8.698e+09 -1.282e+08 -8.698e+04 -1479 8.698e+04 -1282 x7 -2.301e+07 -3.131e+05 2.301e+07 -4.175e+05 -230.1 -3.131 230.1 -4.175 x8 -8.698e+09 -1.282e+08 8.698e+09 -1.479e+08 -8.698e+04 -1282 8.698e+04 -1479 B = u1 u2 u3 u4 x1 0 0 0 0 x2 0 0 0 0 x3 0 0 0 0 x4 0 0 0 0 x5 7384 -1.744e+06 -1846 -8.721e+05 x6 -1.744e+06 5.494e+08 8.721e+05 3.845e+08 x7 -1846 8.721e+05 7384 1.744e+06 x8 -8.721e+05 3.845e+08 1.744e+06 5.494e+08 C = x1 x2 x3 x4 x5 x6 x7 x8 y1 1 0 0 0 0 0 0 0 y2 0 1 0 0 0 0 0 0 y3 0 0 1 0 0 0 0 0 y4 0 0 0 1 0 0 0 0 D = u1 u2 u3 u4 y1 0 0 0 0 y2 0 0 0 0 y3 0 0 0 0 y4 0 0 0 0 Continuous-time state-space model.
step(sys1, 10)
The piezobender is clearly unstable in the open loop.
Calculate the optimal Gain matrix
q = 1.0e0; % to make the response faster, increase exponent of q with r fixed at 1
Q = q*eye(2*n);
r = 1.0e0; % to reduce the energy consumption, increase exponent of r with q fixed at 1
R = r*eye(n);
G = lqr(A, B, Q, R)
G = 4×8
0.5024 -0.0078 0.4976 -0.0078 0.9935 -0.0000 0.0076 -0.0000 -0.1397 0.9415 0.1397 0.0541 -0.0002 1.0000 0.0002 -0.0000 0.4976 0.0078 0.5024 0.0078 0.0076 0.0000 0.9935 0.0000 -0.1397 0.0541 0.1397 0.9415 -0.0002 -0.0000 0.0002 1.0000
Check the step responses of the closed-loop MIMO system
H = A - B*G;
sys2 = ss(H, B, C, D)
sys2 = A = x1 x2 x3 x4 x5 x6 x7 x8 x1 0 0 0 0 1 0 0 0 x2 0 0 0 0 0 1 0 0 x3 0 0 0 0 0 0 1 0 x4 0 0 0 0 0 0 0 1 x5 2.264e+07 2.107e+06 -2.265e+07 1.229e+06 -7634 1.744e+06 2094 8.721e+05 x6 -8.567e+09 -6.86e+08 8.568e+09 -5.2e+08 1.833e+06 -5.493e+08 -9.609e+05 -3.845e+08 x7 -2.265e+07 -1.229e+06 2.264e+07 -2.107e+06 2094 -8.721e+05 -7634 -1.744e+06 x8 -8.568e+09 -5.2e+08 8.567e+09 -6.86e+08 9.609e+05 -3.845e+08 -1.833e+06 -5.493e+08 B = u1 u2 u3 u4 x1 0 0 0 0 x2 0 0 0 0 x3 0 0 0 0 x4 0 0 0 0 x5 7384 -1.744e+06 -1846 -8.721e+05 x6 -1.744e+06 5.494e+08 8.721e+05 3.845e+08 x7 -1846 8.721e+05 7384 1.744e+06 x8 -8.721e+05 3.845e+08 1.744e+06 5.494e+08 C = x1 x2 x3 x4 x5 x6 x7 x8 y1 1 0 0 0 0 0 0 0 y2 0 1 0 0 0 0 0 0 y3 0 0 1 0 0 0 0 0 y4 0 0 0 1 0 0 0 0 D = u1 u2 u3 u4 y1 0 0 0 0 y2 0 0 0 0 y3 0 0 0 0 y4 0 0 0 0 Continuous-time state-space model.
step(sys2, 10)
The piezobender is now stable in closed-loop.
If we define row vectors in the Gain matrix as
,
then we know that
.
.
.
Hence, solving for ,
.
Similarly, we can also get
.
In order to tune the Gain matrix until the desired closed-loop response is obtained, adjust the values of q and r. See comments above in the MATLAB code.
##### 9 个评论显示 7更早的评论隐藏 7更早的评论
Fajwa Noor 2022-6-23

Sorry for opening this topic again, but can i ask if this matlab code implied finite element method?
Hope to hear from you soon. Thank you.
Sam Chak 2022-6-23
It's alright. However, I didn't use any finite element method in the code. I merely used the a kind of matrix algebra manipulation method to obtain the values in the gain matrix G to stabilize the piezobender system as described by the mathematics that you showed to me.
G = lqr(A, B, Q, R)
I also took a look into your recent question:
It seems like the piezobender system is only a part of a much larger system that is modeled in Simscape. The piezobender is connected to the Vibration system on the left and the Full wave rectifier on the right. But I don't see Ctr, Cro, Rtr, Rro in the piezobender dynamics above. Would you clarify? You can also check with someone in your lab who is familiar with the system.

### 类别

Help CenterFile Exchange 中查找有关 Classical Control Design 的更多信息

R2021b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by