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

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.

[Sam Chak]

Can you provide the mathematical model of the energy harvester? It maybe easier to analyze and design a set mathematical law that satisfies the desired optimization criteria and efficiency level. Since there are different operating conditions, then a Gain Scheduling technique or a Variable Structure System are better-suited.

[Fajwa Noor]

https://www.mathworks.com/help/physmod/sps/ref/piezobender.html#:~:text=A%20piezo%20bender%20is%20a%20piezoelectric%20device%20that,material%20with%20a%20polarization%20perpendicular%20to%20the%20stack.

Here in this link it provide the mathematical model for piezo bender.

[Sam Chak]

Thanks. Found it. Now, what type of "sophistication" that you want to design for ?

[Fajwa Noor]

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

Answer 1

[Sam Chak]

Hi @Fajwa Noor

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:

https://www.mathworks.com/help/control/ref/lqr.html

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)

https://www.mathworks.com/help/control/ug/designing-a-family-of-pid-controllers-for-multiple-operating-points.html

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.

[Fajwa Noor]

I see thank you for the suggestion and sharing, actually I'm still quite new and have very little of knowledge of Simulink and Matlab. here, I attach my simulink file that I have add a close loop system in my simulink model using PID Controller for you to review. But the results is nothing much change

[Sam Chak]

Don't mention it, @Fajwa Noor. It is just helping the community. I didn't open the Simulink file yet. If you perform due diligence to extract info from the Simulink file and put the values into the parameters in the matrix differential equation (not a state-space yet),

then I maybe able to show you how to apply the LQR command.

[Fajwa Noor]

where the parameters I extract from the simulink model are:

[Sam Chak]

Hi @Fajwa Noor

Since you already have the numbers, please type out the matrices in MATLAB format.

% 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;

Matrices

m  = (rho*l*w*d)/420

m = 5.1594e-06

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]

M1 = 4×4

156.0000 0.6985 54.0000 -0.4128 0.6985 0.0040 0.4128 -0.0030 54.0000 0.4128 156.0000 -0.6985 -0.4128 -0.0030 -0.6985 0.0040

M  = m*M1

M = 4×4

1.0e-03 * 0.8049 0.0036 0.2786 -0.0021 0.0036 0.0000 0.0021 -0.0000 0.2786 0.0021 0.8049 -0.0036 -0.0021 -0.0000 -0.0036 0.0000

You can check that there are some discrepancies between my M1 and yours.

Your turn to complete the matrice K, B, Kk, and Bb.

To be honest, the dynamics of the piezo bender something new to me. The mass matrix size is . But I can try if it possible to design for .

[Fajwa Noor]

% 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;

m  = (rho*l*w*d)/420

m = 5.1594e-06

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]

M1 = 4×4

156.0000 0.6985 54.0000 -0.4128 0.6985 0.0040 0.4128 -0.0030 54.0000 0.4128 156.0000 -0.6985 -0.4128 -0.0030 -0.6985 0.0040

M  = m*M1

M = 4×4

1.0e-03 * 0.8049 0.0036 0.2786 -0.0021 0.0036 0.0000 0.0021 -0.0000 0.2786 0.0021 0.8049 -0.0036 -0.0021 -0.0000 -0.0036 0.0000

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]

K = 4×4

712.3989 11.3093 -712.3989 11.3093 11.3093 0.2394 -11.3093 0.1197 -712.3989 -11.3093 712.3989 -11.3093 11.3093 0.1197 -11.3093 0.2394

B = kb*K + bm*M 

B = 4×4

0.0071 0.0001 -0.0071 0.0001 0.0001 0.0000 -0.0001 0.0000 -0.0071 -0.0001 0.0071 -0.0001 0.0001 0.0000 -0.0001 0.0000

I only know how to find K and B. How do I find Kk and Bb, please advice.

I see, actually Im doing this for my FYP which is to study the behavior of piezo bender effect in energy harvesting application such as the output power produce at difference frequency.

[Sam Chak]

Hi @Fajwa Noor

I have edited my Answer to include the LQR design. There maybe other ways of designing the closed-loop system effectively and efficiently. Please consult your FYP supervisor on this matter. You can show the LQR design results to your supervisor.

I think this will be the closure of the topic. If you have new questions about implementing LQR in Simulink, please open/ask a new Question and consult the Simulink experts. If you like the design materials and think they are helpful, please Vote the Answer.

[Fajwa Noor]

Hi @Sam Chak

Thank you so much for explaining the closed-loop design as detail as possible, I really appreciate it.

have a nice day ahead. Okay noted i will open a new question

[Fajwa Noor]

Hi @Sam Chak.

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]

Hi @Fajwa Noor

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:

https://www.mathworks.com/matlabcentral/answers/1742510-how-do-i-implement-lqr-linear-quadratic-regulator-design-script-in-my-piezo-bender-simulink-model

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.