From 063d1c79bb1da9821a60dd6afea068e13fbc4a17 Mon Sep 17 00:00:00 2001 From: "americo.cunhajr@gmail.com" <36556019+americocunhajr@users.noreply.github.com> Date: Wed, 17 Jul 2024 16:07:02 -0300 Subject: [PATCH] Update optimization module --- STONEHENGE-1.0/.DS_Store | Bin 10244 -> 10244 bytes .../CrossEntropyOpt_Example.m | 39 ++ .../HarvesterOpt_CrossEntropy2var.m | 419 +++++++++++++++++ .../HarvesterOpt_CrossEntropy2varNoise.m | 419 +++++++++++++++++ .../HarvesterOpt_CrossEntropy4var.m | 399 +++++++++++++++++ .../HarvesterOpt_DirectSearch2var.m | 398 +++++++++++++++++ .../HarvesterOpt_DirectSearch2varNoise.m | 397 +++++++++++++++++ .../HarvesterOpt_DirectSearch4var.m | 335 ++++++++++++++ .../HarvesterOpt-1.0/PiezoMagBeam_ObjFunc.m | 62 +++ .../PiezoMagBeam_ObjFuncNoise.m | 65 +++ .../HarvesterOpt-1.0/PiezoMagBeam_RHS.m | 46 ++ .../graph_binarymap.m | 16 +- .../graph_contour_pnt.m | 20 +- .../graph_type2.m} | 84 ++-- .../plot_ce_animation.m | 14 +- STONEHENGE-1.0/HarvesterOpt-1.0/test01chaos.m | 181 ++++++++ .../trandn.m | 0 .../NonconvexOptimization/.DS_Store | Bin 6148 -> 0 bytes .../main_cross_entropy_example.m | 28 -- .../main_piezomagbeam_opt_ce.m | 403 ----------------- .../main_piezomagbeam_opt_ds.m | 421 ------------------ .../piezomagbeam_opt_PerfFunc.m | 63 --- .../plot_ce_animation_video.m | 121 ----- 23 files changed, 2831 insertions(+), 1099 deletions(-) create mode 100644 STONEHENGE-1.0/HarvesterOpt-1.0/CrossEntropyOpt_Example.m create mode 100644 STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_CrossEntropy2var.m create mode 100644 STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_CrossEntropy2varNoise.m create mode 100644 STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_CrossEntropy4var.m create mode 100644 STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_DirectSearch2var.m create mode 100644 STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_DirectSearch2varNoise.m create mode 100644 STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_DirectSearch4var.m create mode 100644 STONEHENGE-1.0/HarvesterOpt-1.0/PiezoMagBeam_ObjFunc.m create mode 100644 STONEHENGE-1.0/HarvesterOpt-1.0/PiezoMagBeam_ObjFuncNoise.m create mode 100644 STONEHENGE-1.0/HarvesterOpt-1.0/PiezoMagBeam_RHS.m rename STONEHENGE-1.0/{NonconvexOptimization => HarvesterOpt-1.0}/graph_binarymap.m (92%) rename STONEHENGE-1.0/{NonconvexOptimization => HarvesterOpt-1.0}/graph_contour_pnt.m (87%) rename STONEHENGE-1.0/{NonconvexOptimization/graph_contourf_pnt.m => HarvesterOpt-1.0/graph_type2.m} (59%) rename STONEHENGE-1.0/{NonconvexOptimization => HarvesterOpt-1.0}/plot_ce_animation.m (92%) create mode 100644 STONEHENGE-1.0/HarvesterOpt-1.0/test01chaos.m rename STONEHENGE-1.0/{NonconvexOptimization => HarvesterOpt-1.0}/trandn.m (100%) delete mode 100644 STONEHENGE-1.0/NonconvexOptimization/.DS_Store delete mode 100644 STONEHENGE-1.0/NonconvexOptimization/main_cross_entropy_example.m delete mode 100644 STONEHENGE-1.0/NonconvexOptimization/main_piezomagbeam_opt_ce.m delete mode 100644 STONEHENGE-1.0/NonconvexOptimization/main_piezomagbeam_opt_ds.m delete mode 100644 STONEHENGE-1.0/NonconvexOptimization/piezomagbeam_opt_PerfFunc.m delete mode 100644 STONEHENGE-1.0/NonconvexOptimization/plot_ce_animation_video.m diff --git a/STONEHENGE-1.0/.DS_Store b/STONEHENGE-1.0/.DS_Store index 45a0230c4f39ba3c1fbb2b83f71bcd20864eb016..255c6b92efaab64628f57f117e85936d9e50d748 100644 GIT binary patch delta 39 tcmZn(XbG6$F8U^o9{djZwW1_JD?o6W^8a4u|M+{~_E0OFiw0ssZQ3|{~M delta 550 zcmZn(XbG6$&uFkQU^hRb!DI&kmCXhM?5y>o41Nsx40#O6V6u!Mm7#*cpP_)Egdvk5 z7l^AE62WYsnw)gQ;N<+=0tPUUU;$Dj=+4b|aY@R_PXg-VSk}BaBD49BBLhP{X-1@A gF@h8m3NnzrpuoVe*;Mor=Vo>VgVB&M6+O=g09m eps +t = t+1; +x = mu + sigma*randn(N,1); +SX = S(x); +sortSX = sortrows([x SX],2); +Xel = sortSX((N - Nel + 1):N,1); +mu = mean(Xel); +sigma = std(Xel); +fprintf('%g %6.9f %6.9f %6.9f \n', t, S(mu),mu, sigma) +end + +mu +toc diff --git a/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_CrossEntropy2var.m b/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_CrossEntropy2var.m new file mode 100644 index 0000000..ef530d5 --- /dev/null +++ b/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_CrossEntropy2var.m @@ -0,0 +1,419 @@ +% ----------------------------------------------------------------- +% HarvesterOpt_CrossEntropy2var.m +% ----------------------------------------------------------------- +% This program solves an optimization problem that seaks to +% maximize the mean output power of a piezo-magneto-elastic +% energy harvester system, defined by +% +% -- t0 + T +% | +% power := 1/T | lambda v(t)^2 dt +% | +% -- t = t0 +% +% in which the voltage come from the following dynamical system +% +% d2x/dt2 + 2*ksi*dx/dt - 0.5*x*(1-x^2) - chi*v = f*cos(Omega*t) +% dv/dt + lambda*v + kappa*dx/dt = 0 +% + +% initial conditions, +% +% where +% +% x(t) - dimensionless displacement of the beam tip +% v(t) - dimensionless voltage across the load resistance +% t - dimensionless time +% ksi - mechanical damping ratio +% chi - dimensionless piezoeletric coupling term (mechanical) +% f - dimensionless excitation amplitude +% Omega - dimensionless excitation frequency +% lambda - dimensionless time constant reciprocal +% kappa - dimensionless piezoeletric coupling term (eletrical) +% +% The solution strategy is based on the cross-entropy method, +% a stochastic metaheuristic that 'transforms' the optimization +% problem into a rare event estimation problem. +% +% The design variables are: (f,omega) +% +% Reference: +% A. Cunha Jr +% Enhancing the performance of a bistable energy harvesting +% device via the cross-entropy method (2020) +% ----------------------------------------------------------------- +% programmer: Americo Cunha +% americo.cunha@uerj.br +% +% last update: October 23, 2020 +% ----------------------------------------------------------------- + + +clc +clear +close all + + +% program header +% ----------------------------------------------------------- +disp(' ---------------------------------------------------') +disp(' Piezo-Magneto-Elastic Harvester Optimization ') +disp(' (cross-entropy method) ') +disp(' ') +disp(' by ') +disp(' Americo Cunha (UERJ) ') +disp(' americo.cunha@uerj.br ') +disp(' ---------------------------------------------------') +% ----------------------------------------------------------- + + + +% simulation information +% ----------------------------------------------------------- +case_name = 'HarvesterOpt_CE_2var'; + +disp(' '); +disp([' Case Name: ',num2str(case_name)]); +disp(' '); +% ----------------------------------------------------------- + + + +% define physical parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining physical parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +ksi = 0.01; % mechanical damping ratio +chi = 0.05; % dimensionless piezoeletric coupling term (mechanical) +lambda = 0.05; % dimensionless time constant reciprocal +kappa = 0.5; % dimensionless piezoeletric coupling term (eletrical) +%f = 0.115; % dimensionless excitation amplitude +%Omega = 0.8; % dimensionless excitation frequency + +x0 = 1.0; % dimensionless initial displacement +xdot0 = 0.0; % dimensionless initial velocity +v0 = 0.0; % dimensionless initial voltage + +toc +% ----------------------------------------------------------- + + +% 0-1 test for chaos parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining model parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% RNG seed (seed is fixed for reprodutibility) +rng_stream = RandStream('mt19937ar','Seed',30081984); +RandStream.setGlobalStream(rng_stream); + +% tolerance for 0-1 test +tol01 = 0.1; + +% parameter c for 0-1 test +cmin = 0.0; +cmax = 2*pi; + +% number of test repeats +Nc = 100; + +% oversampling flag +OSflag = 0; + +toc +% ----------------------------------------------------------- + + +% define ODE solver parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining ODE solver parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% initial time of analysis +t0 = 0.0; + +% final time of analysis +t1 = 2500.0; + +% time interval of analysis +tspan = [t0 t1]; + +% inicial condition vector +IC = [x0; xdot0; v0]; + +% ODE solver optional parameters +%opt = odeset('RelTol',1.0e-9,'AbsTol',1.0e-6); + +toc +% ----------------------------------------------------------- + + +% solve optimization problem via cross-entropy method +%--------------------------------------------------------------- +tic +disp(' '); +disp(' --- optimization via cross-entropy method --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% number of design variables +Ndv = 2; + +% maximum number of iterations +maxiter = 100; + +% cross-entropy tolerance +tolCE = 1.0e-3; + +% smoothing parameter (0 < alpha <= 1) +% -- set alpha = 1 for no smoothing -- +alpha = 0.7; +one_minus_alpha = 1.0 - alpha; + +% dynamic smoothing parameters +% (0.8 <= beta <= 0.99) +% (q is a interger between 5 and 10) +beta = 0.8; +q = 5; + +% penalty parameter +H = 10.0; + +% elite samples percentage +rho = 0.1; + +% number of cross-entropy samples +NCE = 50; + +% elite samples index +Nelite = NCE - ceil(rho*NCE) + 1; + +% initialize level counter +t = 0; + +% preallocate memory for PerfFunc +S = zeros(NCE,maxiter); + +% preallocate memory for power function +power = zeros(NCE,maxiter); + +% preallocate memory for 0-1 classifier +K01 = zeros(NCE,maxiter); + +% preallocate memory for design variables samples +%X = zeros(NCE*maxiter,Ndv); +X1 = zeros(NCE,maxiter); +X2 = zeros(NCE,maxiter); + +% initialize ObjFunc maximum +S_opt = -Inf; + +% initialize power maximum +power_opt = -Inf; + +% initialize K01 maximum +K01_opt = -Inf; + +% initialize maximum point +X_opt = zeros(Ndv,1); + +% limits for the design variables: (f,Omega) +X_min = [0.08; 0.75]; +X_max = [0.10; 0.85]; + +% initialize mean and std dev vectors + mu = X_min + (X_max-X_min).*rand(Ndv,1); +sigma = 5*(X_max-X_min); + +% initialize old mean and std dev vectors + mu0 = mu; +sigma0 = sigma; + +% inf-norm of sigma +sigma_max = norm(sigma,Inf); + +% hyper parameters +hyperparam.tspan = tspan; +hyperparam.IC = IC; +hyperparam.cmin = cmin; +hyperparam.cmax = cmax; +hyperparam.Nc = Nc; +hyperparam.tol01 = tol01; +hyperparam.OSflag = OSflag; +hyperparam.H = H; + +% parameters +param.ksi = ksi; +param.chi = chi; +param.lambda = lambda; +param.kappa = kappa; + +% define data file name +file_name = [case_name,'.dat']; + +% open data file +fileID = fopen(file_name,'w'); + +% define data format +formatSpec1 = '%03d %+0.4f %0.4f %1.4f %.4f %.4f %.4f %.4f \n'; +formatSpec2 = '%+0.4f %0.4f %1.4f %.4f %.4f \n'; + +% print global optimum in the data file +fprintf(fileID,'\n t S_max P_max K01 mu1 mu2 sigma1 sigma2\n'); + +% display extreme values on screen +fprintf('\n t S_max P_max K01 mu1 mu2 sigma1 sigma2\n'); + +while sigma_max > tolCE && t < maxiter + + % update level counter + t = t + 1; + + % limit vectors for truncated Gaussian + supp_p_l = ((X_min - mu)./sigma)*ones(1,NCE); + supp_p_h = ((X_max - mu)./sigma)*ones(1,NCE); + %supp_X1_l = ((X_min(1) - mu(1))./sigma(1))*ones(NCE,1); + %supp_X1_h = ((X_max(1) - mu(1))./sigma(1))*ones(NCE,1); + %supp_X2_l = ((X_min(2) - mu(2))./sigma(2))*ones(NCE,1); + %supp_X2_h = ((X_max(2) - mu(2))./sigma(2))*ones(NCE,1); + + % generate samples from truncated normal distribution + X1(:,t) = mu(1) + trandn(supp_p_l(1,:),supp_p_h(1,:))*sigma(1); + X2(:,t) = mu(2) + trandn(supp_p_l(2,:),supp_p_h(2,:))*sigma(2); + %X1(:,t) = mu(1) + trandn(supp_X1_l,supp_X1_h)*sigma(1); + %X2(:,t) = mu(2) + trandn(supp_X2_l,supp_X2_h)*sigma(2); + + % evaluate objective function at the samples + for n=1:NCE + + % update parameters + param.f = X1(n,t); + param.Omega = X2(n,t); + + % penalized objective function S(x) + [S(n,t),power(n,t),K01(n,t)] = ... + PiezoMagBeam_ObjFunc(param,hyperparam); + end + + % sort objective function evaluations + [S_sort,I_sort] = sort(S(:,t)); + + % update mean value + mu(1) = mean(X1(I_sort(Nelite:end),t)); + mu(2) = mean(X2(I_sort(Nelite:end),t)); + + % update standard deviation + sigma(1) = std(X1(I_sort(Nelite:end),t)); + sigma(2) = std(X2(I_sort(Nelite:end),t)); + + % estimator for objective function maximum + S_max = S_sort(Nelite); + P_max = power(I_sort(Nelite),t); + K01_max = K01(I_sort(Nelite),t); + + % smoothing + mu = alpha*mu + one_minus_alpha*mu0; + + % update dynamics smoothing parameter + beta_t = beta*(1 - (1-1/t)^q); + + % dynamic smoothing + sigma = beta_t*sigma + (1-beta_t)*sigma0; + + % save old mean + mu0 = mu; + + % save old std dev + sigma0 = sigma; + + % inf-norm of sigma + sigma_max = norm(sigma,Inf); + + % update global maximum + if S_max > S_opt + S_opt = S_max; + power_opt = P_max; + K01_opt = K01_max; + X_opt(1) = X1(I_sort(Nelite),t); + X_opt(2) = X2(I_sort(Nelite),t); + end + + % print local extreme values in a data file + fprintf(fileID,formatSpec1,... + t,S_max,P_max,K01_max,mu(1),mu(2),sigma(1),sigma(2)); + + % print local extreme values on screen + fprintf(formatSpec1,... + t,S_max,P_max,K01_max,mu(1),mu(2),sigma(1),sigma(2)); +end + +% print global optimum in the data file +fprintf(fileID,'\n\nGlobal maximum (ObjFunc power K01 f Omega):\n'); +fprintf(fileID,formatSpec2,... + S_opt,power_opt,K01_opt,X_opt(1),X_opt(2)); +fprintf('\n'); + +% close data file +fclose(fileID); + + +% display global optimum on the screen +fprintf('\n\nGlobal maximum (ObjFunc power K01 f Omega):\n'); +fprintf(formatSpec2,... + S_opt,power_opt,K01_opt,X_opt(1),X_opt(2)); + + +disp(' '); +time_elapsed = toc +%--------------------------------------------------------------- + + +% save simulation results +% ----------------------------------------------------------- +tic +disp(' ') +disp(' --- saving simulation results --- '); +disp(' '); +disp(' ... '); +disp(' '); + +save([num2str(case_name),'.mat']); + +toc +% ----------------------------------------------------------- + + + + +% animate cross-entropy iteration +% ........................................................... +vtitle = 'cross-entropy method'; +xlab = 'excitation amplitude'; +ylab = 'excitation frequency'; +xmin = X_min(1); +xmax = X_max(1); +ymin = X_min(2); +ymax = X_max(2); +xref = 0.0994; % obtained via direct search with a 256 x 256 grid +yref = 0.7771; % obtained via direct search with a 256 x 256 grid +vname = [num2str(case_name),'__animation']; + +plot_ce_animation(X1,X2,(1:t),xref,yref,... + vtitle,xlab,ylab,xmin,xmax,ymin,ymax); +% ........................................................... + + + diff --git a/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_CrossEntropy2varNoise.m b/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_CrossEntropy2varNoise.m new file mode 100644 index 0000000..31a42a4 --- /dev/null +++ b/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_CrossEntropy2varNoise.m @@ -0,0 +1,419 @@ +% ----------------------------------------------------------------- +% HarvesterOpt_CrossEntropy2varNoise.m +% ----------------------------------------------------------------- +% This program solves an optimization problem that seaks to +% maximize the mean output power of a piezo-magneto-elastic +% energy harvester system, defined by +% +% -- t0 + T +% | +% power := 1/T | lambda v(t)^2 dt +% | +% -- t = t0 +% +% in which the voltage come from the following dynamical system +% +% d2x/dt2 + 2*ksi*dx/dt - 0.5*x*(1-x^2) - chi*v = f*cos(Omega*t) +% dv/dt + lambda*v + kappa*dx/dt = 0 +% + +% initial conditions, +% +% where +% +% x(t) - dimensionless displacement of the beam tip +% v(t) - dimensionless voltage across the load resistance +% t - dimensionless time +% ksi - mechanical damping ratio +% chi - dimensionless piezoeletric coupling term (mechanical) +% f - dimensionless excitation amplitude +% Omega - dimensionless excitation frequency +% lambda - dimensionless time constant reciprocal +% kappa - dimensionless piezoeletric coupling term (eletrical) +% +% The solution strategy is based on the cross-entropy method, +% a stochastic metaheuristic that 'transforms' the optimization +% problem into a rare event estimation problem. +% +% The design variables are: (f,omega) +% +% Reference: +% A. Cunha Jr +% Enhancing the performance of a bistable energy harvesting +% device via the cross-entropy method (2020) +% ----------------------------------------------------------------- +% programmer: Americo Cunha +% americo.cunha@uerj.br +% +% last update: October 23, 2020 +% ----------------------------------------------------------------- + + +clc +clear +close all + + +% program header +% ----------------------------------------------------------- +disp(' ---------------------------------------------------') +disp(' Noisily Piezo-Magneto-Elastic Harvester Optimizer ') +disp(' (cross-entropy method) ') +disp(' ') +disp(' by ') +disp(' Americo Cunha (UERJ) ') +disp(' americo.cunha@uerj.br ') +disp(' ---------------------------------------------------') +% ----------------------------------------------------------- + + + +% simulation information +% ----------------------------------------------------------- +case_name = 'HarvesterOpt_CE_2var_Noise'; + +disp(' '); +disp([' Case Name: ',num2str(case_name)]); +disp(' '); +% ----------------------------------------------------------- + + + +% define physical parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining physical parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +ksi = 0.01; % mechanical damping ratio +chi = 0.05; % dimensionless piezoeletric coupling term (mechanical) +lambda = 0.05; % dimensionless time constant reciprocal +kappa = 0.5; % dimensionless piezoeletric coupling term (eletrical) +%f = 0.115; % dimensionless excitation amplitude +%Omega = 0.8; % dimensionless excitation frequency + +x0 = 1.0; % dimensionless initial displacement +xdot0 = 0.0; % dimensionless initial velocity +v0 = 0.0; % dimensionless initial voltage + +toc +% ----------------------------------------------------------- + + +% 0-1 test for chaos parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining model parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% RNG seed (seed is fixed for reprodutibility) +rng_stream = RandStream('mt19937ar','Seed',30081984); +RandStream.setGlobalStream(rng_stream); + +% tolerance for 0-1 test +tol01 = 0.1; + +% parameter c for 0-1 test +cmin = 0.0; +cmax = 2*pi; + +% number of test repeats +Nc = 100; + +% oversampling flag +OSflag = 0; + +toc +% ----------------------------------------------------------- + + +% define ODE solver parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining ODE solver parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% initial time of analysis +t0 = 0.0; + +% final time of analysis +t1 = 2500.0; + +% time interval of analysis +tspan = [t0 t1]; + +% inicial condition vector +IC = [x0; xdot0; v0]; + +% ODE solver optional parameters +%opt = odeset('RelTol',1.0e-9,'AbsTol',1.0e-6); + +toc +% ----------------------------------------------------------- + + +% solve optimization problem via cross-entropy method +%--------------------------------------------------------------- +tic +disp(' '); +disp(' --- optimization via cross-entropy method --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% number of design variables +Ndv = 2; + +% maximum number of iterations +maxiter = 200; + +% cross-entropy tolerance +tolCE = 1/512; + +% smoothing parameter (0 < alpha <= 1) +% -- set alpha = 1 for no smoothing -- +alpha = 0.7; +one_minus_alpha = 1.0 - alpha; + +% dynamic smoothing parameters +% (0.8 <= beta <= 0.99) +% (q is a interger between 5 and 10) +beta = 0.8; +q = 5; + +% penalty parameter +H = 10.0; + +% elite samples percentage +rho = 0.1; + +% number of cross-entropy samples +NCE = 50; + +% elite samples index +Nelite = NCE - ceil(rho*NCE) + 1; + +% initialize level counter +t = 0; + +% preallocate memory for PerfFunc +S = zeros(NCE,maxiter); + +% preallocate memory for power function +power = zeros(NCE,maxiter); + +% preallocate memory for 0-1 classifier +K01 = zeros(NCE,maxiter); + +% preallocate memory for design variables samples +%X = zeros(NCE*maxiter,Ndv); +X1 = zeros(NCE,maxiter); +X2 = zeros(NCE,maxiter); + +% initialize ObjFunc maximum +S_opt = -Inf; + +% initialize power maximum +power_opt = -Inf; + +% initialize K01 maximum +K01_opt = -Inf; + +% initialize maximum point +X_opt = zeros(Ndv,1); + +% limits for the design variables: (f,Omega) +X_min = [0.08; 0.75]; +X_max = [0.10; 0.85]; + +% initialize mean and std dev vectors + mu = X_min + (X_max-X_min).*rand(Ndv,1); +sigma = 5*(X_max-X_min); + +% initialize old mean and std dev vectors + mu0 = mu; +sigma0 = sigma; + +% inf-norm of sigma +sigma_max = norm(sigma,Inf); + +% hyper parameters +hyperparam.tspan = tspan; +hyperparam.IC = IC; +hyperparam.cmin = cmin; +hyperparam.cmax = cmax; +hyperparam.Nc = Nc; +hyperparam.tol01 = tol01; +hyperparam.OSflag = OSflag; +hyperparam.H = H; + +% parameters +param.ksi = ksi; +param.chi = chi; +param.lambda = lambda; +param.kappa = kappa; + +% define data file name +file_name = [case_name,'.dat']; + +% open data file +fileID = fopen(file_name,'w'); + +% define data format +formatSpec1 = '%03d %+0.4f %0.4f %1.4f %.4f %.4f %.4f %.4f \n'; +formatSpec2 = '%+0.4f %0.4f %1.4f %.4f %.4f \n'; + +% print global optimum in the data file +fprintf(fileID,'\n t S_max P_max K01 mu1 mu2 sigma1 sigma2\n'); + +% display extreme values on screen +fprintf('\n t S_max P_max K01 mu1 mu2 sigma1 sigma2\n'); + +while sigma_max > tolCE && t < maxiter + + % update level counter + t = t + 1; + + % limit vectors for truncated Gaussian + supp_p_l = ((X_min - mu)./sigma)*ones(1,NCE); + supp_p_h = ((X_max - mu)./sigma)*ones(1,NCE); + %supp_X1_l = ((X_min(1) - mu(1))./sigma(1))*ones(NCE,1); + %supp_X1_h = ((X_max(1) - mu(1))./sigma(1))*ones(NCE,1); + %supp_X2_l = ((X_min(2) - mu(2))./sigma(2))*ones(NCE,1); + %supp_X2_h = ((X_max(2) - mu(2))./sigma(2))*ones(NCE,1); + + % generate samples from truncated normal distribution + X1(:,t) = mu(1) + trandn(supp_p_l(1,:),supp_p_h(1,:))*sigma(1); + X2(:,t) = mu(2) + trandn(supp_p_l(2,:),supp_p_h(2,:))*sigma(2); + %X1(:,t) = mu(1) + trandn(supp_X1_l,supp_X1_h)*sigma(1); + %X2(:,t) = mu(2) + trandn(supp_X2_l,supp_X2_h)*sigma(2); + + % evaluate objective function at the samples + for n=1:NCE + + % update parameters + param.f = X1(n,t); + param.Omega = X2(n,t); + + % penalized objective function S(x) + [S(n,t),power(n,t),K01(n,t)] = ... + PiezoMagBeam_ObjFuncNoise(param,hyperparam); + end + + % sort objective function evaluations + [S_sort,I_sort] = sort(S(:,t)); + + % update mean value + mu(1) = mean(X1(I_sort(Nelite:end),t)); + mu(2) = mean(X2(I_sort(Nelite:end),t)); + + % update standard deviation + sigma(1) = std(X1(I_sort(Nelite:end),t)); + sigma(2) = std(X2(I_sort(Nelite:end),t)); + + % estimator for objective function maximum + S_max = S_sort(Nelite); + P_max = power(I_sort(Nelite),t); + K01_max = K01(I_sort(Nelite),t); + + % smoothing + mu = alpha*mu + one_minus_alpha*mu0; + + % update dynamics smoothing parameter + beta_t = beta*(1 - (1-1/t)^q); + + % dynamic smoothing + sigma = beta_t*sigma + (1-beta_t)*sigma0; + + % save old mean + mu0 = mu; + + % save old std dev + sigma0 = sigma; + + % inf-norm of sigma + sigma_max = norm(sigma,Inf); + + % update global maximum + if S_max > S_opt + S_opt = S_max; + power_opt = P_max; + K01_opt = K01_max; + X_opt(1) = X1(I_sort(Nelite),t); + X_opt(2) = X2(I_sort(Nelite),t); + end + + % print local extreme values in a data file + fprintf(fileID,formatSpec1,... + t,S_max,P_max,K01_max,mu(1),mu(2),sigma(1),sigma(2)); + + % print local extreme values on screen + fprintf(formatSpec1,... + t,S_max,P_max,K01_max,mu(1),mu(2),sigma(1),sigma(2)); +end + +% print global optimum in the data file +fprintf(fileID,'\n\nGlobal maximum (ObjFunc power K01 f Omega):\n'); +fprintf(fileID,formatSpec2,... + S_opt,power_opt,K01_opt,X_opt(1),X_opt(2)); +fprintf('\n'); + +% close data file +fclose(fileID); + + +% display global optimum on the screen +fprintf('\n\nGlobal maximum (ObjFunc power K01 f Omega):\n'); +fprintf(formatSpec2,... + S_opt,power_opt,K01_opt,X_opt(1),X_opt(2)); + + +disp(' '); +time_elapsed = toc +%--------------------------------------------------------------- + + +% save simulation results +% ----------------------------------------------------------- +tic +disp(' ') +disp(' --- saving simulation results --- '); +disp(' '); +disp(' ... '); +disp(' '); + +save([num2str(case_name),'.mat']); + +toc +% ----------------------------------------------------------- + + + + +% animate cross-entropy iteration +% ........................................................... +vtitle = 'cross-entropy method'; +xlab = 'excitation amplitude'; +ylab = 'excitation frequency'; +xmin = X_min(1); +xmax = X_max(1); +ymin = X_min(2); +ymax = X_max(2); +xref = 0.0998; % obtained via direct search with a 256 x 256 grid +yref = 0.7763; % obtained via direct search with a 256 x 256 grid +vname = [num2str(case_name),'__animation']; + +plot_ce_animation(X1,X2,(1:t),xref,yref,... + vtitle,xlab,ylab,xmin,xmax,ymin,ymax); +% ........................................................... + + + diff --git a/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_CrossEntropy4var.m b/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_CrossEntropy4var.m new file mode 100644 index 0000000..ff344f6 --- /dev/null +++ b/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_CrossEntropy4var.m @@ -0,0 +1,399 @@ +% ----------------------------------------------------------------- +% HarvesterOpt_CrossEntropy4var.m +% ----------------------------------------------------------------- +% This program solves an optimization problem that seaks to +% maximize the mean output power of a piezo-magneto-elastic +% energy harvester system, defined by +% +% -- t0 + T +% | +% power := 1/T | lambda v(t)^2 dt +% | +% -- t = t0 +% +% in which the voltage come from the following dynamical system +% +% d2x/dt2 + 2*ksi*dx/dt - 0.5*x*(1-x^2) - chi*v = f*cos(Omega*t) +% dv/dt + lambda*v + kappa*dx/dt = 0 +% + +% initial conditions, +% +% where +% +% x(t) - dimensionless displacement of the beam tip +% v(t) - dimensionless voltage across the load resistance +% t - dimensionless time +% ksi - mechanical damping ratio +% chi - dimensionless piezoeletric coupling term (mechanical) +% f - dimensionless excitation amplitude +% Omega - dimensionless excitation frequency +% lambda - dimensionless time constant reciprocal +% kappa - dimensionless piezoeletric coupling term (eletrical) +% +% The solution strategy is based on the cross-entropy method, +% a stochastic metaheuristic that 'transforms' the optimization +% problem into a rare event estimation problem. +% +% The design variables are: (ksi,chi,lambda,kappa) +% +% Reference: +% A. Cunha Jr +% Enhancing the performance of a bistable energy harvesting +% device via the cross-entropy method (2020) +% ----------------------------------------------------------------- +% programmer: Americo Cunha +% americo.cunha@uerj.br +% +% last update: October 22, 2020 +% ----------------------------------------------------------------- + + +clc +clear +close all + + +% program header +% ----------------------------------------------------------- +disp(' ---------------------------------------------------') +disp(' Piezo-Magneto-Elastic Harvester Optimization ') +disp(' (cross-entropy method) ') +disp(' ') +disp(' by ') +disp(' Americo Cunha (UERJ) ') +disp(' americo.cunha@uerj.br ') +disp(' ---------------------------------------------------') +% ----------------------------------------------------------- + + + +% simulation information +% ----------------------------------------------------------- +case_name = 'HarvesterOpt_CE_4var'; + +disp(' '); +disp([' Case Name: ',num2str(case_name)]); +disp(' '); +% ----------------------------------------------------------- + + + +% define physical parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining physical parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +%ksi = 0.01; % mechanical damping ratio +%chi = 0.05; % dimensionless piezoeletric coupling term (mechanical) +%lambda = 0.05; % dimensionless time constant reciprocal +%kappa = 0.5; % dimensionless piezoeletric coupling term (eletrical) +f = 0.083; % dimensionless excitation amplitude +Omega = 0.8; % dimensionless excitation frequency + +x0 = 1.0; % dimensionless initial displacement +xdot0 = 0.0; % dimensionless initial velocity +v0 = 0.0; % dimensionless initial voltage + +toc +% ----------------------------------------------------------- + + +% 0-1 test for chaos parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining model parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% tolerance for 0-1 test +tol01 = 0.1; + +% parameter c for 0-1 test +cmin = 0.0; +cmax = 2*pi; + +% number of test repeats +Nc = 100; + +% oversampling flag +OSflag = 0; + +toc +% ----------------------------------------------------------- + + +% define ODE solver parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining ODE solver parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% initial time of analysis +t0 = 0.0; + +% final time of analysis +t1 = 2500.0; + +% time interval of analysis +tspan = [t0 t1]; + +% inicial condition vector +IC = [x0; xdot0; v0]; + +% ODE solver optional parameters +%opt = odeset('RelTol',1.0e-9,'AbsTol',1.0e-6); + +toc +% ----------------------------------------------------------- + + +% solve optimization problem via cross-entropy method +%--------------------------------------------------------------- +tic +disp(' '); +disp(' --- optimization via cross-entropy method --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% number of design variables +Ndv = 4; + +% RNG seed (seed is fixed for reprodutibility) +rng_stream = RandStream('mt19937ar','Seed',30081984); +RandStream.setGlobalStream(rng_stream); + +% maximum number of iterations +maxiter = 100; + +% cross-entropy tolerance +tolCE = 1.0e-2; + +% smoothing parameter (0 < alpha <= 1) +% -- set alpha = 1 for no smoothing -- +alpha = 0.7; +one_minus_alpha = 1.0 - alpha; + +% dynamic smoothing parameters +% (0.8 <= beta <= 0.99) +% (q is a interger between 5 and 10) +beta = 0.8; +q = 5; + +% penalty parameter +H = 10.0; + +% elite samples percentage +rho = 0.1; + +% number of cross-entropy samples +NCE = 100; + +% elite samples index +Nelite = NCE - ceil(rho*NCE) + 1; + +% initialize level counter +t = 0; + +% preallocate memory for PerfFunc +S = zeros(NCE,maxiter); + +% preallocate memory for power function +power = zeros(NCE,maxiter); + +% preallocate memory for 0-1 classifier +K01 = zeros(NCE,maxiter); + +% preallocate memory for design variables samples +X1 = zeros(NCE,maxiter); +X2 = zeros(NCE,maxiter); +X3 = zeros(NCE,maxiter); +X4 = zeros(NCE,maxiter); + +% initialize PerfFunc maximum value +S_opt = -Inf; + +% initialize power maximum value +power_opt = -Inf; + +% initialize K01 maximum value +K01_opt = -Inf; + +% initialize maximum point +X_opt = zeros(Ndv,1); + +% limits for the design variables: (ksi,chi,lambda,kappa) +X_min = [0.01; 0.050; 0.050; 0.50]; +X_max = [0.05; 0.200; 0.200; 1.50]; + +% initialize mean and std dev vectors + mu = X_min + (X_max-X_min).*rand(Ndv,1); +sigma = 5*(X_max-X_min); + +% initialize old mean and std dev vectors + mu0 = mu; +sigma0 = sigma; + +% inf-norm of sigma +sigma_max = norm(sigma,Inf); + +% hyper parameters +hyperparam.tspan = tspan; +hyperparam.IC = IC; +hyperparam.cmin = cmin; +hyperparam.cmax = cmax; +hyperparam.Nc = Nc; +hyperparam.tol01 = tol01; +hyperparam.OSflag = OSflag; +hyperparam.H = H; + +% parameters +param.f = f; +param.Omega = Omega; + +% define data file name +file_name = [case_name,'.dat']; + +% open data file +fileID = fopen(file_name,'w'); + +% define data format +formatSpec1 = '%03d %+0.4f %0.4f %1.4f %.4f %.4f %.4f %.4f %.4f \n'; +formatSpec2 = '%+0.4f %0.4f %1.4f %.4f %.4f %.4f %.4f \n'; + +% print global optimum in the data file +fprintf(fileID,'\n t S_max P_max K01 mu1 mu2 mu3 mu4 sigma_max \n'); + +% display extreme values on screen +fprintf('\n t S_max P_max K01 mu1 mu2 mu3 mu4 sigma_max \n'); + +while sigma_max > tolCE && t < maxiter + + % update level counter + t = t + 1; + + % limit vectors for truncated Gaussian + supp_p_l = ((X_min - mu)./sigma)*ones(1,NCE); + supp_p_h = ((X_max - mu)./sigma)*ones(1,NCE); + + % generate samples from truncated normal distribution + X1(:,t) = mu(1) + trandn(supp_p_l(1,:),supp_p_h(1,:))*sigma(1); + X2(:,t) = mu(2) + trandn(supp_p_l(2,:),supp_p_h(2,:))*sigma(2); + X3(:,t) = mu(3) + trandn(supp_p_l(3,:),supp_p_h(3,:))*sigma(3); + X4(:,t) = mu(4) + trandn(supp_p_l(4,:),supp_p_h(4,:))*sigma(4); + + % evaluate performance function at the samples + for n=1:NCE + + % update parameters + param.ksi = X1(n,t); + param.chi = X2(n,t); + param.lambda = X3(n,t); + param.kappa = X4(n,t); + + % penalized performance function S(x) + [S(n,t),power(n,t),K01(n,t)] = ... + PiezoMagBeam_ObjFunc(param,hyperparam); + end + + % sort performance function evaluations + [S_sort,I_sort] = sort(S(:,t)); + + % update mean value + mu(1) = mean(X1(I_sort(Nelite:end),t)); + mu(2) = mean(X2(I_sort(Nelite:end),t)); + mu(3) = mean(X3(I_sort(Nelite:end),t)); + mu(4) = mean(X4(I_sort(Nelite:end),t)); + + % update standard deviation + sigma(1) = std(X1(I_sort(Nelite:end),t)); + sigma(2) = std(X2(I_sort(Nelite:end),t)); + sigma(3) = std(X3(I_sort(Nelite:end),t)); + sigma(4) = std(X4(I_sort(Nelite:end),t)); + + % estimator for PerfFunc maximum + S_max = S_sort(Nelite); + P_max = power(I_sort(Nelite),t); + K01_max = K01(I_sort(Nelite),t); + + % smoothing + mu = alpha*mu + one_minus_alpha*mu0; + + % update dynamics smoothing parameter + beta_t = beta*(1 - (1-1/t)^q); + + % dynamic smoothing + sigma = beta_t*sigma + (1-beta_t)*sigma0; + + % save old mean + mu0 = mu; + + % save old std dev + sigma0 = sigma; + + % inf-norm of sigma + sigma_max = norm(sigma,Inf); + + % update global maximum + if S_max > S_opt + S_opt = S_max; + power_opt = P_max; + K01_opt = K01_max; + X_opt(1) = X1(I_sort(Nelite),t); + X_opt(2) = X2(I_sort(Nelite),t); + X_opt(3) = X3(I_sort(Nelite),t); + X_opt(4) = X4(I_sort(Nelite),t); + end + + % print local extreme values in a data file + fprintf(fileID,formatSpec1,... + t,S_max,P_max,K01_max,mu(1),mu(2),mu(3),mu(4),sigma_max); + + % print local extreme values on screen + fprintf(formatSpec1,... + t,S_max,P_max,K01_max,mu(1),mu(2),mu(3),mu(4),sigma_max); +end + +% print global optimum in the data file +fprintf(fileID,'\n\nGlobal maximum (ObjFunc power K01 ksi chi lambda kappa):\n'); +fprintf(fileID,formatSpec2,... + S_opt,power_opt,K01_opt,X_opt(1),X_opt(2),X_opt(3),X_opt(4)); +fprintf('\n'); + +% close data file +fclose(fileID); + + +% display global optimum on the screen +fprintf('\n\nGlobal maximum (ObjFunc power K01 ksi chi lambda kappa):\n'); +fprintf(formatSpec2,... + S_opt,power_opt,K01_opt,X_opt(1),X_opt(2),X_opt(3),X_opt(4)); + +disp(' '); +time_elapsed = toc +%--------------------------------------------------------------- + + +% save simulation results +% ----------------------------------------------------------- +tic +disp(' ') +disp(' --- saving simulation results --- '); +disp(' '); +disp(' ... '); +disp(' '); + +save([num2str(case_name),'.mat']); + +toc +% ----------------------------------------------------------- + diff --git a/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_DirectSearch2var.m b/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_DirectSearch2var.m new file mode 100644 index 0000000..3db6ae8 --- /dev/null +++ b/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_DirectSearch2var.m @@ -0,0 +1,398 @@ +% ----------------------------------------------------------------- +% HarvesterOpt_DirectSearch2var.m +% ----------------------------------------------------------------- +% This program solves an optimization problem that seaks to +% maximize the mean output power of a piezo-magneto-elastic +% energy harvester system, defined by +% +% -- t0 + T +% | +% power := 1/T | lambda v(t)^2 dt +% | +% -- t = t0 +% +% in which the voltage come from the following dynamical system +% +% d2x/dt2 + 2*ksi*dx/dt - 0.5*x*(1-x^2) - chi*v = f*cos(Omega*t) +% dv/dt + lambda*v + kappa*dx/dt = 0 +% + +% initial conditions, +% +% where +% +% x(t) - dimensionless displacement of the beam tip +% v(t) - dimensionless voltage across the load resistance +% t - dimensionless time +% ksi - mechanical damping ratio +% chi - dimensionless piezoeletric coupling term (mechanical) +% f - dimensionless excitation amplitude +% Omega - dimensionless excitation frequency +% lambda - dimensionless time constant reciprocal +% kappa - dimensionless piezoeletric coupling term (eletrical) +% +% The solution strategy is based on a direct search in a fine mesh, +% where the optimum is found by exhaustion of the design points. +% +% The design variables are: (f,Omega) +% +% Reference: +% A. Cunha Jr +% Enhancing the performance of a bistable energy harvesting +% device via the cross-entropy method (2020) +% ----------------------------------------------------------------- +% programmer: Americo Cunha +% americo.cunha@uerj.br +% +% last update: October 4, 2020 +% ----------------------------------------------------------------- + +clc +clear +close all + + +% program header +% ----------------------------------------------------------- +disp(' ---------------------------------------------------') +disp(' Piezo-Magneto-Elastic Harvester Optimizator ') +disp(' (direct search) ') +disp(' ') +disp(' by ') +disp(' Americo Cunha (UERJ) ') +disp(' americo.cunha@uerj.br ') +disp(' ---------------------------------------------------') +% ----------------------------------------------------------- + + + +% simulation information +% ----------------------------------------------------------- +case_name = 'HarvesterOpt_DS_2var'; + +disp(' '); +disp([' Case Name: ',num2str(case_name)]); +disp(' '); +% ----------------------------------------------------------- + + + +% define physical parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining model parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +ksi = 0.01; % mechanical damping ratio +chi = 0.05; % dimensionless piezoeletric coupling term (mechanical) +lambda = 0.05; % dimensionless time constant reciprocal +kappa = 0.5; % dimensionless piezoeletric coupling term (eletrical) +%f = 0.115; % dimensionless excitation amplitude +%Omega = 0.8; % dimensionless excitation frequency + +x0 = 1.0; % dimensionless initial displacement +xdot0 = 0.0; % dimensionless initial velocity +v0 = 0.0; % dimensionless initial voltage + +toc +%--------------------------------------------------------------- + + +% 0-1 test for chaos parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining 0-1 test for chaos parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% RNG seed (seed is fixed for reprodutibility) +rng_stream = RandStream('mt19937ar','Seed',30081984); +RandStream.setGlobalStream(rng_stream); + +% tolerance for 0-1 test +tol01 = 0.1; + +% parameter c for 0-1 test +cmin = 0.0; +cmax = 2*pi; + +% number of test repeats +Nc = 100; + +% oversampling flag +OSflag = 0; + +toc +% ----------------------------------------------------------- + + +% define ODE solver parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining ODE solver parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% initial time of analysis +t0 = 0.0; + +% final time of analysis +t1 = 2500.0; + +% time interval of analysis +tspan = [t0 t1]; + +% inicial condition vector +IC = [x0; xdot0; v0]; + +% ODE solver optional parameters +%opt = odeset('RelTol',1.0e-9,'AbsTol',1.0e-6); + +toc +% ----------------------------------------------------------- + + + +% solve optimization problem via direct search +%--------------------------------------------------------------- +tic +disp(' '); +disp(' --- optimization via direct search --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% number of design variables +Ndv = 2; + +% number of points for design variables discretization +Np1 = 256; +Np2 = 256; + +% number of points in the numerical grid +Ngrid = Np1*Np2; + + +% square root of Ngrid +sqrt_Ngrid = round(sqrt(Ngrid)); + +% penalty parameter +H = 10.0; + +% limits for the design variables: (f,Omega) +X_min = [0.08 0.75]; +X_max = [0.10 0.85]; + +% numerical grid for domain discretization +X1 = linspace(X_min(1),X_max(1),Np1); +X2 = linspace(X_min(2),X_max(2),Np2); + +% preallocate memory for the ObjFunc +S = zeros(Ngrid,1); + +% preallocate memory for power function +% (x = columns and y = lines) +power = zeros(Ngrid,1); + +% preallocate memory for 0-1 classifier +% (x = columns and y = lines) +K01 = zeros(Ngrid,1); + +% preallocate memory for the optimum point +X_opt = zeros(Ndv,1); + +% initialize ObjFunc maximum +S_opt = -Inf; + +% initialize power maximum +power_opt = -Inf; + +% initialize K01 maximum +K01_opt = -Inf; + +% hyper parameters +hyperparam.tspan = tspan; +hyperparam.IC = IC; +hyperparam.cmin = cmin; +hyperparam.cmax = cmax; +hyperparam.Nc = Nc; +hyperparam.tol01 = tol01; +hyperparam.OSflag = OSflag; +hyperparam.H = H; + +% parameters +param.ksi = ksi; +param.chi = chi; +param.lambda = lambda; +param.kappa = kappa; + +% define data file name +file_name = [case_name,'.dat']; + +% open data file +fileID = fopen(file_name,'w'); + +% define data format +formatSpec = '%+0.4f %0.4f %1.4f %.4f %.4f \n'; + +% print global optimum in the data file +fprintf(fileID,'\n S_max P_max K01 mu1 mu2 \n'); + +for np1 = 1:Np1 + + % update design variable + f = X1(np1); + + for np2 = 1:Np2 + + ngrid = np2 + (np1-1)*Np2; + + % print loop indicador + if mod(ngrid,sqrt_Ngrid) == 0 + disp(['ngrid = ',num2str(ngrid)]); + end + + % update design variable + Omega = X2(np2); + + % define physical paramters vector + param.f = f; + param.Omega = Omega; + + % penalized objective function S(x) + [S(ngrid),power(ngrid),K01(ngrid)] = ... + PiezoMagBeam_ObjFunc(param,hyperparam); + + % update global maximum + if S(ngrid) > S_opt + + S_opt = S(ngrid); + power_opt = power(ngrid); + K01_opt = K01(ngrid); + X_opt = [f; Omega]; + end + + % print local values in data file + fprintf(fileID,formatSpec,... + S(ngrid),power(ngrid),K01(ngrid),X1(np1),X2(np2)); + + end +end + +% print global optimum in the data file +%fprintf(fileID,'\n'); +fprintf(fileID,'\n\nGlobal maximum (ObjFunc power K01 f Omega):\n'); +fprintf(fileID,formatSpec,S_opt,power_opt,K01_opt,X_opt(1),X_opt(2)); + +% close data file +fclose(fileID); + +% display global optimum on the screen +fprintf('\n\nGlobal maximum (ObjFunc power K01 f Omega):\n'); +fprintf(formatSpec,S_opt,power_opt,K01_opt,X_opt(1),X_opt(2)); +fprintf('\n'); + +disp(' '); +time_elapsed = toc +%--------------------------------------------------------------- + + + +% save simulation results +% ----------------------------------------------------------- +tic +disp(' ') +disp(' --- saving simulation results --- '); +disp(' '); +disp(' ... '); +disp(' '); + +save([num2str(case_name),'.mat']); + +toc +% ----------------------------------------------------------- + + +% post processing +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- post processing --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% map to a 2d data structure +% (x = columns and y = lines) + S = reshape(S ,Np2,Np1); +power = reshape(power,Np2,Np1); + K01 = reshape(K01 ,Np2,Np1); + + +% plot performance function countourmap +% ...................................................... +gtitle = ' objective function contour map'; +xlab = ' excitation amplitude'; +ylab = ' excitation frequency'; +xmin = X_min(1); +xmax = X_max(1); +ymin = X_min(2); +ymax = X_max(2); +gname = [case_name,'__objfunc']; +flag = 'eps'; +fig1 = graph_contour_pnt(X1,X2,S'/S_opt,... + X_opt(1),X_opt(2),... + gtitle,xlab,ylab,xmin,xmax,ymin,ymax,gname,flag); +%close(fig1); +% ...................................................... + + +% plot mean power countourmap +% ...................................................... +gtitle = ' mean power contour map'; +xlab = ' excitation amplitude'; +ylab = ' excitation frequency'; +xmin = X_min(1); +xmax = X_max(1); +ymin = X_min(2); +ymax = X_max(2); +gname = [case_name,'__mean_power']; +flag = 'eps'; +fig2 = graph_contour_pnt(X1,X2,power'/power_opt,... + X_opt(1),X_opt(2),... + gtitle,xlab,ylab,xmin,xmax,ymin,ymax,gname,flag); + +%close(fig2); +% ...................................................... + + +% plot 0-1 identification parameter +% ...................................................... +xlab = 'excitation amplitude'; +ylab = 'excitation frequency'; +label0 = 'regular'; +label1 = 'chaos'; +gtitle = '0-1 test for chaos'; +xmin = X_min(1); +xmax = X_max(1); +ymin = X_min(2); +ymax = X_max(2); +zmin = 0; +zmax = 1; +gname = [case_name,'__01test']; +flag = 'eps'; +fig3 = graph_binarymap(X1,X2,K01,... + X_opt(1),X_opt(2),... + gtitle,xlab,ylab,label0,label1,... + xmin,xmax,ymin,ymax,zmin,zmax,gname,flag); +%close(fig3); +% ...................................................... + +toc +% ----------------------------------------------------------- diff --git a/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_DirectSearch2varNoise.m b/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_DirectSearch2varNoise.m new file mode 100644 index 0000000..6dd2c1f --- /dev/null +++ b/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_DirectSearch2varNoise.m @@ -0,0 +1,397 @@ +% ----------------------------------------------------------------- +% HarvesterOpt_DirectSearch2varNoise.m +% ----------------------------------------------------------------- +% This program solves an optimization problem that seaks to +% maximize the mean output power of a piezo-magneto-elastic +% energy harvester system, defined by +% +% -- t0 + T +% | +% power := 1/T | lambda v(t)^2 dt +% | +% -- t = t0 +% +% in which the voltage come from the following dynamical system +% +% d2x/dt2 + 2*ksi*dx/dt - 0.5*x*(1-x^2) - chi*v = f*cos(Omega*t) +% dv/dt + lambda*v + kappa*dx/dt = 0 +% + +% initial conditions, +% +% where +% +% x(t) - dimensionless displacement of the beam tip +% v(t) - dimensionless voltage across the load resistance +% t - dimensionless time +% ksi - mechanical damping ratio +% chi - dimensionless piezoeletric coupling term (mechanical) +% f - dimensionless excitation amplitude +% Omega - dimensionless excitation frequency +% lambda - dimensionless time constant reciprocal +% kappa - dimensionless piezoeletric coupling term (eletrical) +% +% The solution strategy is based on a direct search in a fine mesh, +% where the optimum is found by exhaustion of the design points. +% +% The design variables are: (f,Omega) +% +% Reference: +% A. Cunha Jr +% Enhancing the performance of a bistable energy harvesting +% device via the cross-entropy method (2020) +% ----------------------------------------------------------------- +% programmer: Americo Cunha +% americo.cunha@uerj.br +% +% last update: October 23, 2020 +% ----------------------------------------------------------------- + +clc +clear +close all + + +% program header +% ----------------------------------------------------------- +disp(' ---------------------------------------------------') +disp(' Noisily Piezo-Magneto-Elastic Harvester Optimizer ') +disp(' (direct search) ') +disp(' ') +disp(' by ') +disp(' Americo Cunha (UERJ) ') +disp(' americo.cunha@uerj.br ') +disp(' ---------------------------------------------------') +% ----------------------------------------------------------- + + + +% simulation information +% ----------------------------------------------------------- +case_name = 'HarvesterOpt_DS_2var_Noise'; + +disp(' '); +disp([' Case Name: ',num2str(case_name)]); +disp(' '); +% ----------------------------------------------------------- + + + +% define physical parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining model parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +ksi = 0.01; % mechanical damping ratio +chi = 0.05; % dimensionless piezoeletric coupling term (mechanical) +lambda = 0.05; % dimensionless time constant reciprocal +kappa = 0.5; % dimensionless piezoeletric coupling term (eletrical) +%f = 0.115; % dimensionless excitation amplitude +%Omega = 0.8; % dimensionless excitation frequency + +x0 = 1.0; % dimensionless initial displacement +xdot0 = 0.0; % dimensionless initial velocity +v0 = 0.0; % dimensionless initial voltage + +toc +%--------------------------------------------------------------- + + +% 0-1 test for chaos parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining 0-1 test for chaos parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% RNG seed (seed is fixed for reprodutibility) +rng_stream = RandStream('mt19937ar','Seed',30081984); +RandStream.setGlobalStream(rng_stream); + +% tolerance for 0-1 test +tol01 = 0.1; + +% parameter c for 0-1 test +cmin = 0.0; +cmax = 2*pi; + +% number of test repeats +Nc = 100; + +% oversampling flag +OSflag = 0; + +toc +% ----------------------------------------------------------- + + +% define ODE solver parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining ODE solver parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% initial time of analysis +t0 = 0.0; + +% final time of analysis +t1 = 2500.0; + +% time interval of analysis +tspan = [t0 t1]; + +% inicial condition vector +IC = [x0; xdot0; v0]; + +% ODE solver optional parameters +%opt = odeset('RelTol',1.0e-9,'AbsTol',1.0e-6); + +toc +% ----------------------------------------------------------- + + + +% solve optimization problem via direct search +%--------------------------------------------------------------- +tic +disp(' '); +disp(' --- optimization via direct search --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% number of design variables +Ndv = 2; + +% number of points for design variables discretization +Np1 = 256; +Np2 = 256; + +% number of points in the numerical grid +Ngrid = Np1*Np2; + +% square root of Ngrid +sqrt_Ngrid = round(sqrt(Ngrid)); + +% penalty parameter +H = 10.0; + +% limits for the design variables: (f,Omega) +X_min = [0.08 0.75]; +X_max = [0.10 0.85]; + +% numerical grid for domain discretization +X1 = linspace(X_min(1),X_max(1),Np1); +X2 = linspace(X_min(2),X_max(2),Np2); + +% preallocate memory for the ObjFunc +S = zeros(Ngrid,1); + +% preallocate memory for power function +% (x = columns and y = lines) +power = zeros(Ngrid,1); + +% preallocate memory for 0-1 classifier +% (x = columns and y = lines) +K01 = zeros(Ngrid,1); + +% preallocate memory for the optimum point +X_opt = zeros(Ndv,1); + +% initialize ObjFunc maximum +S_opt = -Inf; + +% initialize power maximum +power_opt = -Inf; + +% initialize K01 maximum +K01_opt = -Inf; + +% hyper parameters +hyperparam.tspan = tspan; +hyperparam.IC = IC; +hyperparam.cmin = cmin; +hyperparam.cmax = cmax; +hyperparam.Nc = Nc; +hyperparam.tol01 = tol01; +hyperparam.OSflag = OSflag; +hyperparam.H = H; + +% parameters +param.ksi = ksi; +param.chi = chi; +param.lambda = lambda; +param.kappa = kappa; + +% define data file name +file_name = [case_name,'.dat']; + +% open data file +fileID = fopen(file_name,'w'); + +% define data format +formatSpec = '%+0.4f %0.4f %1.4f %.4f %.4f \n'; + +% print global optimum in the data file +fprintf(fileID,'\n S_max P_max K01 mu1 mu2 \n'); + +for np1 = 1:Np1 + + % update design variable + f = X1(np1); + + for np2 = 1:Np2 + + ngrid = np2 + (np1-1)*Np2; + + % print loop indicador + if mod(ngrid,sqrt_Ngrid) == 0 + disp(['ngrid = ',num2str(ngrid)]); + end + + % update design variable + Omega = X2(np2); + + % define physical paramters vector + param.f = f; + param.Omega = Omega; + + % penalized objective function S(x) + [S(ngrid),power(ngrid),K01(ngrid)] = ... + PiezoMagBeam_ObjFuncNoise(param,hyperparam); + + % update global maximum + if S(ngrid) > S_opt + + S_opt = S(ngrid); + power_opt = power(ngrid); + K01_opt = K01(ngrid); + X_opt = [f; Omega]; + end + + % print local values in data file + fprintf(fileID,formatSpec,... + S(ngrid),power(ngrid),K01(ngrid),X1(np1),X2(np2)); + + end +end + +% print global optimum in the data file +%fprintf(fileID,'\n'); +fprintf(fileID,'\n\nGlobal maximum (ObjFunc power K01 f Omega):\n'); +fprintf(fileID,formatSpec,S_opt,power_opt,K01_opt,X_opt(1),X_opt(2)); + +% close data file +fclose(fileID); + +% display global optimum on the screen +fprintf('\n\nGlobal maximum (ObjFunc power K01 f Omega):\n'); +fprintf(formatSpec,S_opt,power_opt,K01_opt,X_opt(1),X_opt(2)); +fprintf('\n'); + +disp(' '); +time_elapsed = toc +%--------------------------------------------------------------- + + + +% save simulation results +% ----------------------------------------------------------- +tic +disp(' ') +disp(' --- saving simulation results --- '); +disp(' '); +disp(' ... '); +disp(' '); + +save([num2str(case_name),'.mat']); + +toc +% ----------------------------------------------------------- + + +% post processing +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- post processing --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% map to a 2d data structure +% (x = columns and y = lines) + S = reshape(S ,Np2,Np1); +power = reshape(power,Np2,Np1); + K01 = reshape(K01 ,Np2,Np1); + + +% plot performance function countourmap +% ...................................................... +gtitle = ' objective function contour map'; +xlab = ' excitation amplitude'; +ylab = ' excitation frequency'; +xmin = X_min(1); +xmax = X_max(1); +ymin = X_min(2); +ymax = X_max(2); +gname = [case_name,'__objfunc']; +flag = 'eps'; +fig1 = graph_contour_pnt(X1,X2,S'/S_opt,... + X_opt(1),X_opt(2),... + gtitle,xlab,ylab,xmin,xmax,ymin,ymax,gname,flag); +%close(fig1); +% ...................................................... + + +% plot mean power countourmap +% ...................................................... +gtitle = ' mean power contour map'; +xlab = ' excitation amplitude'; +ylab = ' excitation frequency'; +xmin = X_min(1); +xmax = X_max(1); +ymin = X_min(2); +ymax = X_max(2); +gname = [case_name,'__mean_power']; +flag = 'eps'; +fig2 = graph_contour_pnt(X1,X2,power'/power_opt,... + X_opt(1),X_opt(2),... + gtitle,xlab,ylab,xmin,xmax,ymin,ymax,gname,flag); + +%close(fig2); +% ...................................................... + + +% plot 0-1 identification parameter +% ...................................................... +xlab = 'excitation amplitude'; +ylab = 'excitation frequency'; +label0 = 'regular'; +label1 = 'chaos'; +gtitle = '0-1 test for chaos'; +xmin = X_min(1); +xmax = X_max(1); +ymin = X_min(2); +ymax = X_max(2); +zmin = 0; +zmax = 1; +gname = [case_name,'__01test']; +flag = 'eps'; +fig3 = graph_binarymap(X1,X2,K01,... + X_opt(1),X_opt(2),... + gtitle,xlab,ylab,label0,label1,... + xmin,xmax,ymin,ymax,zmin,zmax,gname,flag); +%close(fig3); +% ...................................................... + +toc +% ----------------------------------------------------------- diff --git a/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_DirectSearch4var.m b/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_DirectSearch4var.m new file mode 100644 index 0000000..a2fb94c --- /dev/null +++ b/STONEHENGE-1.0/HarvesterOpt-1.0/HarvesterOpt_DirectSearch4var.m @@ -0,0 +1,335 @@ +% ----------------------------------------------------------------- +% HarvesterOpt_DirectSearch4var.m +% ----------------------------------------------------------------- +% This program solves an optimization problem that seaks to +% maximize the mean output power of a piezo-magneto-elastic +% energy harvester system, defined by +% +% -- t0 + T +% | +% power := 1/T | lambda v(t)^2 dt +% | +% -- t = t0 +% +% in which the voltage come from the following dynamical system +% +% d2x/dt2 + 2*ksi*dx/dt - 0.5*x*(1-x^2) - chi*v = f*cos(Omega*t) +% dv/dt + lambda*v + kappa*dx/dt = 0 +% + +% initial conditions, +% +% where +% +% x(t) - dimensionless displacement of the beam tip +% v(t) - dimensionless voltage across the load resistance +% t - dimensionless time +% ksi - mechanical damping ratio +% chi - dimensionless piezoeletric coupling term (mechanical) +% f - dimensionless excitation amplitude +% Omega - dimensionless excitation frequency +% lambda - dimensionless time constant reciprocal +% kappa - dimensionless piezoeletric coupling term (eletrical) +% +% The solution strategy is based on a direct search in a fine mesh, +% where the optimum is found by exhaustion of the design points. +% +% The design variables are: (ksi,chi,lambda,kappa) +% +% Reference: +% A. Cunha Jr +% Enhancing the performance of a bistable energy harvesting +% device via the cross-entropy method (2020) +% ----------------------------------------------------------------- +% programmer: Americo Cunha +% americo.cunha@uerj.br +% +% last update: October 23, 2020 +% ----------------------------------------------------------------- + +clc +clear +close all + + +% program header +% ----------------------------------------------------------- +disp(' ---------------------------------------------------') +disp(' Piezo-Magneto-Elastic Harvester Optimizer ') +disp(' (direct search) ') +disp(' ') +disp(' by ') +disp(' Americo Cunha (UERJ) ') +disp(' americo.cunha@uerj.br ') +disp(' ---------------------------------------------------') +% ----------------------------------------------------------- + + + +% simulation information +% ----------------------------------------------------------- +case_name = 'HarvesterOpt_DS_4var'; + +disp(' '); +disp([' Case Name: ',num2str(case_name)]); +disp(' '); +% ----------------------------------------------------------- + + + +% define physical parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining model parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +%ksi = 0.01; % mechanical damping ratio +%chi = 0.05; % dimensionless piezoeletric coupling term (mechanical) +%lambda = 0.05; % dimensionless time constant reciprocal +%kappa = 0.5; % dimensionless piezoeletric coupling term (eletrical) +f = 0.083; % dimensionless excitation amplitude +Omega = 0.8; % dimensionless excitation frequency + +x0 = 1.0; % dimensionless initial displacement +xdot0 = 0.0; % dimensionless initial velocity +v0 = 0.0; % dimensionless initial voltage + +toc +%--------------------------------------------------------------- + + +% 0-1 test for chaos parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining 0-1 test for chaos parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% RNG seed (seed is fixed for reprodutibility) +rng_stream = RandStream('mt19937ar','Seed',30081984); +RandStream.setGlobalStream(rng_stream); + +% tolerance for 0-1 test +tol01 = 0.1; + +% parameter c for 0-1 test +cmin = 0.0; +cmax = 2*pi; + +% number of test repeats +Nc = 100; + +% oversampling flag +OSflag = 0; + +toc +% ----------------------------------------------------------- + + +% define ODE solver parameters +% ----------------------------------------------------------- +tic +disp(' '); +disp(' --- defining ODE solver parameters --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% initial time of analysis +t0 = 0.0; + +% final time of analysis +t1 = 2500.0; + +% time interval of analysis +tspan = [t0 t1]; + +% inicial condition vector +IC = [x0; xdot0; v0]; + +% ODE solver optional parameters +%opt = odeset('RelTol',1.0e-9,'AbsTol',1.0e-6); + +toc +% ----------------------------------------------------------- + + + +% solve optimization problem via direct search +%--------------------------------------------------------------- +tic +disp(' '); +disp(' --- optimization via direct search --- '); +disp(' '); +disp(' ... '); +disp(' '); + +% number of design variables +Ndv = 4; + +% number of points for design variables discretization +Np1 = 16; +Np2 = 16; +Np3 = 16; +Np4 = 16; + +% number of points in the numerical grid +Ngrid = Np1*Np2*Np3*Np4; + +% square root of Ngrid +sqrt_Ngrid = round(sqrt(Ngrid)); + +% penalty parameter +H = 10.0; + +% limits for the design variables: (ksi,chi,lambda,kappa) +X_min = [0.01; 0.050; 0.050; 0.50]; +X_max = [0.05; 0.200; 0.200; 1.50]; + +% numerical grid for domain discretization +X1 = linspace(X_min(1),X_max(1),Np1); +X2 = linspace(X_min(2),X_max(2),Np2); +X3 = linspace(X_min(3),X_max(3),Np3); +X4 = linspace(X_min(4),X_max(4),Np4); + +% preallocate memory for PerfFunc +S = zeros(Ngrid,1); + +% preallocate memory for power function +power = zeros(Ngrid,1); + +% preallocate memory for 0-1 classifier +K01 = zeros(Ngrid,1); + +% preallocate memory for the optimum point +X_opt = zeros(Ndv,1); + +% initialize ObjFunc maximum +S_opt = -Inf; + +% initialize power maximum +power_opt = -Inf; + +% initialize K01 maximum +K01_opt = -Inf; + +% hyper parameters +hyperparam.tspan = tspan; +hyperparam.IC = IC; +hyperparam.cmin = cmin; +hyperparam.cmax = cmax; +hyperparam.Nc = Nc; +hyperparam.tol01 = tol01; +hyperparam.OSflag = OSflag; +hyperparam.H = H; + +% parameters +param.f = f; +param.Omega = Omega; + +% define data file name +file_name = [case_name,'.dat']; + +% open data file +fileID = fopen(file_name,'w'); + +% define data format +formatSpec = '%+0.4f %0.4f %1.4f %.4f %.4f %.4f %.4f \n'; + +% print global optimum in the data file +fprintf(fileID,'\n S_max P_max K01 mu1 mu2 mu3 mu4 \n'); + +for np1 = 1:Np1 + + % update design variable + ksi = X1(np1); + + for np2 = 1:Np2 + + % update design variable + chi = X2(np2); + + for np3 = 1:Np3 + + % update design variable + lambda = X3(np3); + + for np4 = 1:Np4 + + ngrid = np4 + (np3-1)*Np4 + ... + (np2-1)*Np4*Np3 + ... + (np1-1)*Np4*Np3*Np2; + + % print loop indicador + if mod(ngrid,sqrt_Ngrid) == 0 + disp(['ngrid = ',num2str(ngrid)]); + end + + % update design variable + kappa = X4(np4); + + % define physical paramters vector + param.ksi = ksi; + param.chi = chi; + param.lambda = lambda; + param.kappa = kappa; + + % penalized objective function S(x) + [S(ngrid),power(ngrid),K01(ngrid)] = ... + PiezoMagBeam_ObjFunc(param,hyperparam); + + % update global maximum + if S(ngrid) > S_opt + + S_opt = S(ngrid); + power_opt = power(ngrid); + K01_opt = K01(ngrid); + X_opt = [ksi; chi; lambda; kappa]; + end + + % print local values in data file + fprintf(fileID,formatSpec,... + S(ngrid),power(ngrid),K01(ngrid),... + X1(np1),X2(np2),X3(np3),X4(np4)); + end + end + end +end + +% print global optimum in the data file +fprintf(fileID,'\n\nGlobal maximum (ObjFunc power K01 ksi chi lambda kappa):\n'); +fprintf(fileID,formatSpec,... + S_opt,power_opt,K01_opt,X_opt(1),X_opt(2),X_opt(3),X_opt(4)); + +% close data file +fclose(fileID); + +% display global optimum on the screen +fprintf('\n\nGlobal maximum (ObjFunc power K01 ksi chi lambda kappa):\n'); +fprintf(formatSpec,... + S_opt,power_opt,K01_opt,X_opt(1),X_opt(2),X_opt(3),X_opt(4)); +fprintf('\n'); + +disp(' '); +time_elapsed = toc +%--------------------------------------------------------------- + + + +% save simulation results +% ----------------------------------------------------------- +tic +disp(' ') +disp(' --- saving simulation results --- '); +disp(' '); +disp(' ... '); +disp(' '); + +save([num2str(case_name),'.mat']); + +toc +% ----------------------------------------------------------- diff --git a/STONEHENGE-1.0/HarvesterOpt-1.0/PiezoMagBeam_ObjFunc.m b/STONEHENGE-1.0/HarvesterOpt-1.0/PiezoMagBeam_ObjFunc.m new file mode 100644 index 0000000..7b189fa --- /dev/null +++ b/STONEHENGE-1.0/HarvesterOpt-1.0/PiezoMagBeam_ObjFunc.m @@ -0,0 +1,62 @@ +% ----------------------------------------------------------------- +% PiezoMagBeam_ObjFunc.m +% ----------------------------------------------------------------- +% Objective function for the optimization problem. +% ----------------------------------------------------------------- +% programmer: Americo Cunha +% americo.cunha@uerj.br +% +% last update: October 23, 2020 +% ----------------------------------------------------------------- + +% ----------------------------------------------------------------- +function [S,power,K01] = PiezoMagBeam_ObjFunc(param,hyperparam) + +% hyper parameters +tspan = hyperparam.tspan; +IC = hyperparam.IC; +cmin = hyperparam.cmin; +cmax = hyperparam.cmax; +Nc = hyperparam.Nc; +tol01 = hyperparam.tol01; +OSflag = hyperparam.OSflag; +H = hyperparam.H; + +% parameters +%f = param.f; +%Omega = param.Omega; +%ksi = param.ksi; +%chi = param.chi; +lambda = param.lambda; +%kappa = param.kappa; + + +% integrate the dynamical system with RK-45 +[time,y] = ode45(@(t,y)PiezoMagBeam_RHS(t,y,param),tspan,IC); + +% number of time steps +Ndt = length(time); + +% number of steps for steady state +Nss = round(0.5*Ndt); + +% number of jumps +Njump = round((Ndt-Nss)/(0.05*Ndt)); + +% temporal interval of analysis +T = time(end)-time(Nss); + +% compute the mean output power +power = (lambda/T)*trapz(time(Nss:end),y(Nss:end,3).^2); + +% apply 0-1 test for chaos in voltage time series +K01 = test01chaos(y(Nss:Njump:end,3),cmin,cmax,Nc,OSflag); + +% penalization function +Penalty = max(K01 - tol01,0); + +% penalized objective function +S = power - H*Penalty; + +end +% ----------------------------------------------------------------- diff --git a/STONEHENGE-1.0/HarvesterOpt-1.0/PiezoMagBeam_ObjFuncNoise.m b/STONEHENGE-1.0/HarvesterOpt-1.0/PiezoMagBeam_ObjFuncNoise.m new file mode 100644 index 0000000..e45c593 --- /dev/null +++ b/STONEHENGE-1.0/HarvesterOpt-1.0/PiezoMagBeam_ObjFuncNoise.m @@ -0,0 +1,65 @@ +% ----------------------------------------------------------------- +% PiezoMagBeam_ObjFuncNoise.m +% ----------------------------------------------------------------- +% Objective function for the optimization problem. +% ----------------------------------------------------------------- +% programmer: Americo Cunha +% americo.cunha@uerj.br +% +% last update: October 23, 2020 +% ----------------------------------------------------------------- + +% ----------------------------------------------------------------- +function [S,power,K01] = PiezoMagBeam_ObjFuncNoise(param,hyperparam) + +% hyper parameters +tspan = hyperparam.tspan; +IC = hyperparam.IC; +cmin = hyperparam.cmin; +cmax = hyperparam.cmax; +Nc = hyperparam.Nc; +tol01 = hyperparam.tol01; +OSflag = hyperparam.OSflag; +H = hyperparam.H; + +% parameters +%f = param.f; +%Omega = param.Omega; +%ksi = param.ksi; +%chi = param.chi; +lambda = param.lambda; +%kappa = param.kappa; + + +% integrate the dynamical system with RK-45 +[time,y] = ode45(@(t,y)PiezoMagBeam_RHS(t,y,param),tspan,IC); + +% number of time steps +Ndt = length(time); + +% number of steps for steady state +Nss = round(0.5*Ndt); + +% number of jumps +Njump = round((Ndt-Nss)/(0.05*Ndt)); + +% temporal interval of analysis +T = time(end)-time(Nss); + +% disturb voltage with noise +y(:,3) = y(:,3) + 0.05*max(y(:,3))*randn(Ndt,1); + +% compute the mean output power +power = (lambda/T)*trapz(time(Nss:end),y(Nss:end,3).^2); + +% apply 0-1 test for chaos in voltage time series +K01 = test01chaos(y(Nss:Njump:end,3),cmin,cmax,Nc,OSflag); + +% penalization function +Penalty = max(K01 - tol01,0); + +% penalized objective function +S = power - H*Penalty; + +end +% ----------------------------------------------------------------- diff --git a/STONEHENGE-1.0/HarvesterOpt-1.0/PiezoMagBeam_RHS.m b/STONEHENGE-1.0/HarvesterOpt-1.0/PiezoMagBeam_RHS.m new file mode 100644 index 0000000..eed1c33 --- /dev/null +++ b/STONEHENGE-1.0/HarvesterOpt-1.0/PiezoMagBeam_RHS.m @@ -0,0 +1,46 @@ +% ----------------------------------------------------------------- +% piezomagbeam.m +% ----------------------------------------------------------------- +% This function defines the right hand side of the follwoing +% nonlinear system of ordinary differential equations +% +% dy1/dt = y2 +% dy2/dt = -2*ksi.y2 + 0.5*y1(1-y1^2) + chi.y3 + f.cos(Omega.t) +% dy3/dt = - lambda.y3 - kappa.y2, +% +% that models the dynamics of a piezo-magneto-elastic beam. +% +% Reference: +% A. Erturk, J. Hoffmann, and D. J. Inman +% A piezomagnetoelastic structure for broadband vibration +% energy harvesting +% Applied Physics Letters +% vol. 94 pp. 254102, 2009 +% ----------------------------------------------------------------- +% programmer: Americo Cunha +% americo.cunha@uerj.br +% +% last update: October 23, 2020 +% ----------------------------------------------------------------- + +% ----------------------------------------------------------------- +function ydot = PiezoMagBeam_RHS(t,y,param) + +% preallocate memory for Ydot vector +ydot = zeros(3,1); + +% physical parameters +f = param.f; +Omega = param.Omega; +ksi = param.ksi; +chi = param.chi; +lambda = param.lambda; +kappa = param.kappa; + +% state space system of equations +ydot(1) = y(2); +ydot(2) = -2*ksi*y(2) + 0.5*y(1)*(1.0-y(1)^2) + chi*y(3) + f*cos(Omega*t); +ydot(3) = -lambda*y(3) - kappa*y(2); + +end +% ----------------------------------------------------------------- diff --git a/STONEHENGE-1.0/NonconvexOptimization/graph_binarymap.m b/STONEHENGE-1.0/HarvesterOpt-1.0/graph_binarymap.m similarity index 92% rename from STONEHENGE-1.0/NonconvexOptimization/graph_binarymap.m rename to STONEHENGE-1.0/HarvesterOpt-1.0/graph_binarymap.m index e874950..1d1c43a 100644 --- a/STONEHENGE-1.0/NonconvexOptimization/graph_binarymap.m +++ b/STONEHENGE-1.0/HarvesterOpt-1.0/graph_binarymap.m @@ -1,7 +1,6 @@ - % ----------------------------------------------------------------- % graph_binarymap.m -% +% ----------------------------------------------------------------- % This functions plots the contour map of a binary scalar % function F: R^2 -> {0,1}. % @@ -25,11 +24,11 @@ % % output: % gname.eps - output file in eps format (optional) -% ----------------------------------------------------------------- -% programmer: Americo Barbosa da Cunha Junior -% americo.cunhajr@gmail.com +% ----------------------------------------------------------------- +% programmer: Americo Cunha +% americo.cunha@uerj.br % -% last update: Dec 14, 2016 +% last update: March 31, 2020 % ----------------------------------------------------------------- % ----------------------------------------------------------------- @@ -67,8 +66,8 @@ set(gca,'XColor',[.3 .3 .3],'YColor',[.3 .3 .3]); set(gca,'YDir','normal') set(gca,'FontName','Helvetica'); - set(gca,'FontSize',18); - set(fh2,'LineWidth',2.0); + set(gca,'FontSize',22); + set(fh2,'LineWidth',3.0); set(fh2,'MarkerSize',20.0); %set(gca,'XTick',xmin:xmax); %set(gca,'YTick',ymin:ymax); @@ -93,7 +92,6 @@ if ( strcmp(flag,'eps') ) saveas(gcf,gname,'epsc2'); %gname = [gname, '.eps']; - %graph_fixPSlinestyle(gname,gname); end diff --git a/STONEHENGE-1.0/NonconvexOptimization/graph_contour_pnt.m b/STONEHENGE-1.0/HarvesterOpt-1.0/graph_contour_pnt.m similarity index 87% rename from STONEHENGE-1.0/NonconvexOptimization/graph_contour_pnt.m rename to STONEHENGE-1.0/HarvesterOpt-1.0/graph_contour_pnt.m index 91f5551..bb46909 100644 --- a/STONEHENGE-1.0/NonconvexOptimization/graph_contour_pnt.m +++ b/STONEHENGE-1.0/HarvesterOpt-1.0/graph_contour_pnt.m @@ -1,7 +1,6 @@ - % ----------------------------------------------------------------- % graph_contour_pnt.m -% +% ----------------------------------------------------------------- % This functions plots the contour map of a scalar function % F: R^2 -> R. % @@ -21,11 +20,11 @@ % % output: % gname.eps - output file in eps format (optional) -% ----------------------------------------------------------------- -% programmer: Americo Barbosa da Cunha Junior -% americo.cunhajr@gmail.com +% ----------------------------------------------------------------- +% programmer: Americo Cunha +% americo.cunha@uerj.br % -% last update: Oct 29, 2014 +% last update: March 31, 2020 % ----------------------------------------------------------------- % ----------------------------------------------------------------- @@ -64,8 +63,8 @@ set(gca,'XColor',[.3 .3 .3],'YColor',[.3 .3 .3]); set(gca,'YDir','normal') set(gca,'FontName','Helvetica'); - set(gca,'FontSize',18); - set(fh2,'LineWidth',2.0); + set(gca,'FontSize',22); + set(fh2,'LineWidth',3.0); set(fh2,'MarkerSize',20.0); %set(gca,'XTick',xmin:xmax); %set(gca,'YTick',ymin:ymax); @@ -80,8 +79,8 @@ else ylim([ymin ymax]); end - labX = xlabel(xlab,'FontSize',18,'FontName','Helvetica'); - labY = ylabel(ylab,'FontSize',18,'FontName','Helvetica'); + labX = xlabel(xlab,'FontSize',22,'FontName','Helvetica'); + labY = ylabel(ylab,'FontSize',22,'FontName','Helvetica'); %set(labX,'interpreter','latex'); %set(labY,'interpreter','latex'); @@ -92,7 +91,6 @@ if ( strcmp(flag,'eps') ) saveas(gcf,gname,'epsc2'); %gname = [gname, '.eps']; - %graph_fixPSlinestyle(gname,gname); end return diff --git a/STONEHENGE-1.0/NonconvexOptimization/graph_contourf_pnt.m b/STONEHENGE-1.0/HarvesterOpt-1.0/graph_type2.m similarity index 59% rename from STONEHENGE-1.0/NonconvexOptimization/graph_contourf_pnt.m rename to STONEHENGE-1.0/HarvesterOpt-1.0/graph_type2.m index f4573b1..71f8121 100644 --- a/STONEHENGE-1.0/NonconvexOptimization/graph_contourf_pnt.m +++ b/STONEHENGE-1.0/HarvesterOpt-1.0/graph_type2.m @@ -1,15 +1,17 @@ % ----------------------------------------------------------------- -% graph_contourf_pnt.m +% graph_type2.m % -% This functions plots the contour map of a scalar function -% F: R^2 -> R. +% This functions plots a graph with two curves. % % input: -% x - x mesh vector -% y - y mesh vector -% F - scalar field +% x1 - x data vector 1 +% y1 - y data vector 1 +% x2 - x data vector 2 +% y2 - y data vector 2 % gtitle - graph title +% leg1 - legend 1 +% leg2 - legend 2 % xlab - x axis label % ylab - y axis label % xmin - x axis minimum value @@ -25,71 +27,83 @@ % programmer: Americo Barbosa da Cunha Junior % americo.cunhajr@gmail.com % -% last update: Oct 29, 2014 +% last update: Mar 21, 2014 % ----------------------------------------------------------------- % ----------------------------------------------------------------- -function fig = graph_contourf_pnt(x,y,F,Xmax,Ymax,gtitle,xlab,ylab,... - xmin,xmax,ymin,ymax,gname,flag) - +function fig = graph_type2(x1,y1,x2,y2,... + gtitle,leg1,leg2,... + xlab,ylab,xmin,xmax,ymin,ymax,gname,flag) + % check number of arguments - if nargin < 13 + if nargin < 14 error('Too few inputs.') - elseif nargin > 14 + elseif nargin > 15 error('Too many inputs.') - elseif nargin == 13 + elseif nargin == 14 flag = 'none'; end + + % check arguments + if length(x1) ~= length(y1) + error('x1 and y1 vectors must be same length') + end + + if length(x2) ~= length(y2) + error('x2 and y2 vectors must be same length') + end + fig = figure('Name',gname,'NumberTitle','off'); - % generate 2D mesh grid - [Xq,Yq] = meshgrid(x,y); - fh = contourf(Xq,Yq,F); - %colormap jet; - %colormap cool - %colormap pink - %colormap hot - colormap bone - colorbar - %scaxis([min(min(F)),max(max(F))]) - hold on - fh2 = plot(Xmax,Ymax,'xb'); - hold off + fh1 = plot(x1,y1,'-b'); + hold all + fh2 = plot(x2,y2,'-r'); set(gcf,'color','white'); - set(gca,'position',[0.15 0.2 0.7 0.7]); + set(gca,'position',[0.2 0.2 0.7 0.7]); set(gca,'Box','on'); set(gca,'TickDir','out','TickLength',[.02 .02]); set(gca,'XMinorTick','on','YMinorTick','on'); - set(gca,'XGrid','off','YGrid','off'); + set(gca,'XGrid','off','YGrid','on'); set(gca,'XColor',[.3 .3 .3],'YColor',[.3 .3 .3]); - set(gca,'YDir','normal') set(gca,'FontName','Helvetica'); set(gca,'FontSize',18); - set(fh2,'LineWidth',2.0); - set(fh2,'MarkerSize',10.0); %set(gca,'XTick',xmin:xmax); %set(gca,'YTick',ymin:ymax); %axis([xmin xmax ymin ymax]); + if ( strcmp(xmin,'auto') || strcmp(xmax,'auto') ) xlim('auto'); else xlim([xmin xmax]); end + if ( strcmp(ymin,'auto') || strcmp(ymax,'auto') ) ylim('auto'); else ylim([ymin ymax]); end - labX = xlabel(xlab,'FontSize',18,'FontName','Helvetica'); - labY = ylabel(ylab,'FontSize',18,'FontName','Helvetica'); + + set(fh1,'LineWidth',1.0); + set(fh1,'MarkerSize',4.0); + set(fh1,'MarkerFaceColor','w'); + set(fh1,'MarkerEdgeColor','b'); + set(fh2,'LineWidth',2.0); + set(fh2,'MarkerSize',4.0); + set(fh2,'MarkerFaceColor','w'); + set(fh2,'MarkerEdgeColor','r'); + leg = legend(leg1,leg2); + set(leg,'FontSize',12); + %set(leg,'interpreter', 'latex'); + labX = xlabel(xlab,'FontSize',16,'FontName','Helvetica'); + labY = ylabel(ylab,'FontSize',16,'FontName','Helvetica'); %set(labX,'interpreter','latex'); %set(labY,'interpreter','latex'); - title(gtitle,'FontSize',20,'FontName','Helvetica'); + hold off + title(gtitle,'FontSize',20,'FontName','Helvetica'); - if ( strcmp(flag,'eps') ) saveas(gcf,gname,'epsc2'); %gname = [gname, '.eps']; diff --git a/STONEHENGE-1.0/NonconvexOptimization/plot_ce_animation.m b/STONEHENGE-1.0/HarvesterOpt-1.0/plot_ce_animation.m similarity index 92% rename from STONEHENGE-1.0/NonconvexOptimization/plot_ce_animation.m rename to STONEHENGE-1.0/HarvesterOpt-1.0/plot_ce_animation.m index 9f2c1ae..a5f2cca 100644 --- a/STONEHENGE-1.0/NonconvexOptimization/plot_ce_animation.m +++ b/STONEHENGE-1.0/HarvesterOpt-1.0/plot_ce_animation.m @@ -1,7 +1,6 @@ - % ----------------------------------------------------------------- % plot_ce_animation.m -% +% ----------------------------------------------------------------- % This functions plots an animation of cross-entropy method % for optimization % @@ -18,11 +17,11 @@ % xmax - x axis maximum value % ymin - y axis minimum value % ymax - y axis maximum value -% ----------------------------------------------------------------- -% programmer: Americo Barbosa da Cunha Junior -% americo.cunhajr@gmail.com +% ----------------------------------------------------------------- +% programmer: Americo Cunha +% americo.cunha@uerj.br % -% last update: Feb 10, 2016 +% last update: October 4, 2020 % ----------------------------------------------------------------- % ----------------------------------------------------------------- @@ -96,10 +95,9 @@ % define video title title(vtitle,'FontSize',16,'FontName','Helvetica'); - gname = ['piezomagbeam_opt_ce_penal_10_NCE_50_maxiter_100__fig',num2str(n)]; + gname = ['animation_',num2str(n)]; saveas(gcf,gname,'epsc2'); gname = [gname, '.eps']; - graph_fixPSlinestyle(gname,gname); pause(1.5) %pause diff --git a/STONEHENGE-1.0/HarvesterOpt-1.0/test01chaos.m b/STONEHENGE-1.0/HarvesterOpt-1.0/test01chaos.m new file mode 100644 index 0000000..9c2c9df --- /dev/null +++ b/STONEHENGE-1.0/HarvesterOpt-1.0/test01chaos.m @@ -0,0 +1,181 @@ +% ----------------------------------------------------------------- +% test01chaos.m +% ----------------------------------------------------------------- +% This functions applies Gottwald-Melbourne 0-1 test for chaos +% to a given time series. Result is near to 0 for non-chaotic +% data and near 1 for chaotic data. +% +% Adapted from Z1TEST.m by Paul Matthews, July 2009 +% http://arxiv.org/pdf/0906.1418v1 +% https://goo.gl/vWSqXs +% +% References: +% Georg A. Gottwald and Ian Melbourne +% The 0-1 Test for Chaos: A review +% In Chaos Detection and Predictability, +% Springer Lecture Notes in Physics 915 +% Editors: C. Skokos, G.A. Gottwald and J. Laskar, 2016 +% http://www.springer.com/gp/book/9783662484081 +% +% D. Bernardini and G. Litak +% An overview of 0-1 test for chaos +% J Braz. Soc. Mech. Sci. Eng. (2016) 38: 1433. +% doi:10.1007/s40430-015-0453-y +% +% input: +% x - (N x 1) observable time series +% cmin - lower bound for test parameter c +% cmax - upper bound for test parameter c +% Nc - number of test repeats +% plotflag - diagnostic plots flag (1 - turn on / 0 - turn off) +% OSflag - oversampling msg flag (1 - turn on / 0 - turn off) +% +% output: +% kmedian - ( 1 x 1) classifier (correlation vector mediam) +% kcorr - (Nc x 1) correlation vector +% c - (Nc x 1) parameter c vector +% p - ( N x Nc) generalized dispalcement p +% q - ( N x Nc) generalized dispalcement q +% OS_ID1 - oversampling indicator 1 +% OS_ID2 - oversampling indicator 2 +% ----------------------------------------------------------------- +% programmer: Americo Cunha +% americo.cunha@uerj.br +% +% last update: October 23, 2020 +% ----------------------------------------------------------------- + +% ----------------------------------------------------------------- +function [kmedian,kcorr,c,p,q,OS_ID1,OS_ID2] = ... + test01chaos(x,cmin,cmax,Nc,OSflag,plotflag) + +% check number of arguments +if nargin < 1 + error('Too few inputs.') +elseif nargin > 6 + error('Too many inputs.') +elseif nargin == 1 + cmin = 0.0; + cmax = 2*pi; + Nc = 100; + OSflag = 1; + plotflag = 0; +elseif nargin == 2 + error('[cmin,cmax] interval should be prescribed.') +elseif nargin == 3 + Nc = 100; + OSflag = 1; + plotflag = 0; +elseif nargin == 4 + OSflag = 1; + plotflag = 0; +elseif nargin == 5 + plotflag = 0; +end + +% check if c interval is valid +if cmax <= cmin || cmin < 0.0 || cmax > 2*pi + error('[cmin,cmax] must be a subset of [0,2*pi] interval.') +end + +% compute the time series dimenions +s = size(x); + +% convert time series to a column vector (if necessary) +if s(1) == 1 + x = x'; +end + +% number of points in time series +N = length(x); + +% number of points in short time series +N0 = round(N/10); + +% check time series size +if N0 < 1 + %error('time series is too short') + N0 = N; +end + +% index vetor for original time series +% (N x 1) vector +j = (1:N)'; + +% random values for c in [cmin,cmax] \subset [0,2*pi] +% (Nc x 1) vector +c = cmin + rand(Nc,1)*(cmax-cmin); + +% square of observable time average +AVG2 = mean(x)^2; + +% reduced time series to represent the line y = x +% (N0 x 1) vector +bisec_line = (1:N0)'; + +% preallocate memory for generalized coordinates +p = zeros(N,Nc); +q = zeros(N,Nc); + +% preallocate memory for mean square displacement +% (N0 x 1) vector +MSD = zeros(N0,1); + +% prealloctate memory for correlation vector (classifier) +% (Nc x 1) vector +kcorr = zeros(Nc,1); + +% loop over samples of parameter c +for it=1:Nc + + % transform time series from (x,xdot) to (p,q) space + p(:,it) = cumsum(x.*cos(j*c(it,1))); + q(:,it) = cumsum(x.*sin(j*c(it,1))); + + % mean square displacement minus oscilatory part + for n=1:N0; + + + MSD(n,1) = mean((p(n+1:N)-p(1:N-n)).^2 + ... + (q(n+1:N)-q(1:N-n)).^2) - ... + AVG2*(1-cos(n*c(it,1)))/(1-cos(c(it,1))); + end + + % correlation coefficient + kcorr(it,1) = corr(bisec_line,MSD); +end + +% correlation vector mediam (classifier) +kmedian = abs(median(kcorr)); + +% two crude attempts to check for oversampling +OS_ID1 = (max(x)-min(x))/mean(abs(diff(x))); +OS_ID2 = median(kcorr(cmean(c))); + +% oversampling message +if OSflag == 1 + + if OS_ID1 > 10 || OS_ID2 > 0.5 + disp('Warning: data is probably oversampled.') + disp('Use coarser sampling or reduce the maximum value of c.') + end +end + +% useful diagnostic plots +if plotflag == 1 + + figure(1); + plot(c,abs(kcorr),'*'); + xlabel('c');ylabel('Kc'); + + figure(2) + plot(bisec_line,MSD); + xlabel('N');ylabel('MSD'); + + figure(3) + plot(p(:,1),q(:,1),'o-r','MarkerEdgeColor','b'); + xlabel('p');ylabel('q'); +end + +return +% ----------------------------------------------------------------- diff --git a/STONEHENGE-1.0/NonconvexOptimization/trandn.m b/STONEHENGE-1.0/HarvesterOpt-1.0/trandn.m similarity index 100% rename from STONEHENGE-1.0/NonconvexOptimization/trandn.m rename to STONEHENGE-1.0/HarvesterOpt-1.0/trandn.m diff --git a/STONEHENGE-1.0/NonconvexOptimization/.DS_Store b/STONEHENGE-1.0/NonconvexOptimization/.DS_Store deleted file mode 100644 index 406d34d5989ffede239bf4559132756e0f4f5110..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHK!A=4(5Pek+L@y+I@VGyq$RDf(4<u^2|n-NG3NBP>6P7lglv*`oAX&dWyUR9c)^VRJ~C8s!#zJ_e8&Ea{RXS$ z)GDso@3W40J32naK+gz?nlf*KF6}eNGuDjfN6+fiafVA=potS)F=DZQh#Kd0O#evF zoFOZZnAy>JOXZxNjQXnc8E8~>tref!idW=_wfBEu^I=2zTxZSijINwRPF;S30p8gn z<=u#GodIXS8Tep8&WD65m>L!l_2{6{Cje1?BUIRyw}j**hN)o@kt384Q;9J(*b^hf zbovvEOAU*NF&)7kAHlvH?1d8a>#Uz>I6@lHtux>ZEHbbiyA`?rPrrZvFD7}<8E^(R zih)pSv>J6T$?n#r#mQZpQ17TB5?4gj6dJpZwINsW0ab-#i8P3*VG)rYivI{G4Q`x) HKV{$z4abjV diff --git a/STONEHENGE-1.0/NonconvexOptimization/main_cross_entropy_example.m b/STONEHENGE-1.0/NonconvexOptimization/main_cross_entropy_example.m deleted file mode 100644 index 8976135..0000000 --- a/STONEHENGE-1.0/NonconvexOptimization/main_cross_entropy_example.m +++ /dev/null @@ -1,28 +0,0 @@ - -clc -clear all -close all - -S = inline('exp(-(x-2).^2) + 0.8*exp(-(x+2).^2)'); - -mu = -6; -sigma = 10; -Nel = 10; -N = 100; -eps = 1E-8; -tic -t=0; - -while sigma > eps -t = t+1; -x = mu + sigma*randn(N,1); -SX = S(x); -sortSX = sortrows([x SX],2); -Xel = sortSX((N - Nel + 1):N,1); -mu = mean(Xel); -sigma = std(Xel); -fprintf('%g %6.9f %6.9f %6.9f \n', t, S(mu),mu, sigma) -end - -mu -toc \ No newline at end of file diff --git a/STONEHENGE-1.0/NonconvexOptimization/main_piezomagbeam_opt_ce.m b/STONEHENGE-1.0/NonconvexOptimization/main_piezomagbeam_opt_ce.m deleted file mode 100644 index d8bcdf3..0000000 --- a/STONEHENGE-1.0/NonconvexOptimization/main_piezomagbeam_opt_ce.m +++ /dev/null @@ -1,403 +0,0 @@ - -% ----------------------------------------------------------------- -% main_piezomagbeam_opt_ce.m -% -% This is the main file for a program which solve an optimization -% problem, via cross-entropy method, that seaks to maximize -% output power of a piezo-magneto-elastic beam, which evolves -% according to -% -% d2x/dt2 + 2*ksi*dx/dt - 0.5*x*(1-x^2) - chi*v = f*cos(Omega*t) -% -% dv/dt + lambda*v + kappa*dx/dt = 0 -% -% + -% -% initial conditions, -% -% where -% -% x(t) - dimensionless displacement of the beam tip -% v(t) - dimensionless voltage across the load resistance -% t - dimensionless time -% ksi - mechanical damping ratio -% chi - dimensionless piezoeletric coupling term (mechanical) -% f - dimensionless excitation amplitude -% Omega - dimensionless excitation frequency -% lambda - dimensionless time constant reciprocal -% kappa - dimensionless piezoeletric coupling term (eletrical) -% -% Reference: -% -% A. Cunha Jr -% Enhancing the performance of a bistable energy harvesting -% device via the cross-entropy method -% 2020 -% ----------------------------------------------------------------- -% programmer: Americo Barbosa da Cunha Junior -% americo.cunhajr@gmail.com -% -% last update: March 6, 2017 -% ----------------------------------------------------------------- - - -clc -clear all -close all - - -% program header -% ----------------------------------------------------------- -disp(' ') -disp(' Piezo-Magneto-Elastic Beam Dynamics ') -disp(' (optimization via cross-entropy) ') -disp(' ') -disp(' by ') -disp(' Americo Barbosa da Cunha Junior ') -disp(' americo.cunhajr@gmail.com ') -disp(' ') -% ----------------------------------------------------------- - - - -% define physical parameters -% ----------------------------------------------------------- -disp(' '); -disp(' --- defining physical parameters --- '); -disp(' '); -disp(' ... '); -disp(' '); - -ksi = 0.01; % mechanical damping ratio -chi = 0.05; % dimensionless piezoeletric coupling term (mechanical) -lambda = 0.05; % dimensionless time constant reciprocal -kappa = 0.5; % dimensionless piezoeletric coupling term (eletrical) -f = 0.083; % dimensionless excitation amplitude -Omega = 0.08; % dimensionless excitation frequency - -x0 = 1.0; % dimensionless initial displacement -xdot0 = 0.0; % dimensionless initial velocity -v0 = 0.0; % dimensionless initial voltage -% ----------------------------------------------------------- - - -% 0-1 test for chaos parameters -% ----------------------------------------------------------- - -% tolerance for 0-1 test -tol01 = 0.1; - -% parameter c for 0-1 test -cmin = 0.0; -cmax = 2*pi; - -% number of test repeats -Nc = 100; - -% oversampling flag -OSflag = 0; -% ----------------------------------------------------------- - - -% ODE solver parameters -% ----------------------------------------------------------- - -% ODE solver optional parameters -%opt = odeset('RelTol',1.0e-9,'AbsTol',1.0e-6); - -% inicial condition vector -IC = [x0; xdot0; v0]; - -% initial time of analysis -t0 = 0.0; - -% final time of analysis -t1 = 2500.0; - -% time interval of analysis -tspan = [t0 t1]; -% ----------------------------------------------------------- - - -% solve optimization problem via cross-entropy method -%--------------------------------------------------------------- -tic - -disp(' '); -disp(' --- optimization problem via cross-entropy --- '); -disp(' '); -disp(' ... '); -disp(' '); - - -% RNG seed -rng_stream = RandStream('mt19937ar','Seed',30081984); -RandStream.setDefaultStream(rng_stream); % Matlab 2009 -%RandStream.setGlobalStream(rng_stream); % Matlab 2013 - -% maximum number of iterations -maxiter = 100; - -% cross-entropy tolerance -tolCE = 1.0e-3; - -% fixed smoothing parameter -% (0 <= alpha <= 1) typically 0.7 -% alpha = 1 for no smoothing -alpha = 0.7; -one_minus_alpha = 1.0 - alpha; - -% dynamic smoothing parameters -% (0.8 <= beta <= 0.99) -% (q is a interger between 5 and 10) -beta = 0.8; -q = 5; - -% penalty parameter -Hpenalty = 10.0; - -% elite samples percentage -rho = 0.1; - -% number of cross-entropy samples -NCE = 25; - -% elite samples index -Nelite = NCE - ceil(rho*NCE) + 1; -%Nelite = NCE -10+1; - -% initialize level counter -t = 0; - -% preallocate memory for PerfFunc -S = zeros(NCE,maxiter); - -% preallocate memory for power function -power = zeros(NCE,maxiter); - -% preallocate memory for 0-1 classifier -K01 = zeros(NCE,maxiter); - -% preallocate memory for p1 samples (excitation amplitude) -p1 = zeros(NCE,maxiter); - -% preallocate memory for p2 samples (excitation frequency) -p2 = zeros(NCE,maxiter); - -% initialize PerfFunc maximum value -S_max_overall = -Inf; - -% initialize power maximum value -power_max_overall = -Inf; - -% initialize K01 maximum value -K01_max_overall = -Inf; - -% initialize PerfFunc maximum point -p1_max_overall = 0.0; -p2_max_overall = 0.0; - -% control parameter 1 -%p1_min = 0.08; %camplitude f -%p1_max = 0.10; % amplitude f -p1_min = 0.04; % coupling chi -p1_max = 0.06; % coupling chi - -% control parameter 2 -%p2_min = 0.75; % frequency w -%p2_max = 0.85; % frequency w -p2_min = 0.4; % coupling kappa -p2_max = 0.6; % coupling kappa - -% support for initial parameter vector -pmin = [p1_min p2_min]; -pmax = [p1_max p2_max]; - -% initialize mean and std dev vectors - mu = pmin + (pmax-pmin).*rand(1,2); -sigma = 5*(pmax-pmin); - -% initialize old mean and std dev vectors - mu0 = mu; -sigma0 = sigma; - -% simulation information -%case_name = ['piezomagbeam_opt_ce_',... -% 'penal_',num2str(Hpenalty),'_',... -% 'NCE_',num2str(NCE),'_',... -% 'maxiter_',num2str(maxiter)]; - -% simulation information -case_name = ['piezomagbeam_optdesign_ce_',... - 'penal_',num2str(Hpenalty),'_',... - 'NCE_',num2str(NCE),'_',... - 'maxiter_',num2str(maxiter)]; - -% define data file name -file_name = [case_name,'.dat']; - -% open data file -fileID = fopen(file_name,'w'); - -% define data format -formatSpec1 = '%03d %+0.4E %+0.4E %+0.4E %0.4E %0.4E \n'; - -% display extreme values on screen -%fprintf('\n t S_max mu1 mu2 sigma1 sigma2\n'); - -while norm(sigma,Inf) > tolCE && t < maxiter - - % update level counter - t = t+1; - - % generate samples from normal distribution - %p1(:,t) = mu(1,1) + randn(NCE,1)*sigma(1,1); - %p2(:,t) = mu(1,2) + randn(NCE,1)*sigma(1,2); - - % limit vectors for truncated Gaussian - supp_p1_l = ((pmin(1,1) - mu(1,1))./sigma(1,1))*ones(NCE,1); - supp_p1_h = ((pmax(1,1) - mu(1,1))./sigma(1,1))*ones(NCE,1); - supp_p2_l = ((pmin(1,2) - mu(1,2))./sigma(1,2))*ones(NCE,1); - supp_p2_h = ((pmax(1,2) - mu(1,2))./sigma(1,2))*ones(NCE,1); - - % generate samples from truncated normal distribution - p1(:,t) = mu(1,1) + trandn(supp_p1_l,supp_p1_h)*sigma(1,1); - p2(:,t) = mu(1,2) + trandn(supp_p2_l,supp_p2_h)*sigma(1,2); - - % evaluate performance function at the samples - for n=1:NCE - - % update parameters - %f = p1(n,t); - %Omega = p2(n,t); - chi = p1(n,t); - kappa = p2(n,t); - - % define physical paramters vector - phys_param = [ksi chi f Omega lambda kappa]; - - % penalized performance function S(x) - [S(n,t),power(n,t),K01(n,t)] = ... - piezomagbeam_opt_PerfFunc(phys_param,tspan,IC,... - cmin,cmax,Nc,tol01,OSflag,Hpenalty,... - p1_min,p1_max,p2_min,p2_max); - end - - % sort performance function evaluations - [S_sort,I_sort] = sort(S(:,t)); - - % update mean value - mu(1,1) = mean(p1(I_sort(Nelite:end),t)); - mu(1,2) = mean(p2(I_sort(Nelite:end),t)); - - % update standard deviation - sigma(1,1) = std(p1(I_sort(Nelite:end),t)); - sigma(1,2) = std(p2(I_sort(Nelite:end),t)); - - % estimator for PerfFunc maximum - S_max = S_sort(Nelite); - - % fixed smoothing - mu = alpha*mu + one_minus_alpha*mu0; - %sigma = alpha*sigma + one_minus_alpha*sigma0; - - % update dynamics smoothing parameter - beta_t = beta*(1 - (1-1/t)^q); - - % dynamic smoothing - sigma = beta_t*sigma + (1-beta_t)*sigma0; - - % save old mean - mu0 = mu; - - % save old std dev - sigma0 = sigma; - - % update global maximum - if S_max > S_max_overall - S_max_overall = S_max; - power_max_overall = power(I_sort(Nelite),t); - K01_max_overall = K01(I_sort(Nelite),t); - p1_max_overall = p1(I_sort(Nelite),t); - p2_max_overall = p2(I_sort(Nelite),t); - end - - % print local extreme values on screen - fprintf(formatSpec1,... - t,S_max,mu(1),mu(2),sigma(1),sigma(2)); - - % print local extreme values in a data file - fprintf(fileID,formatSpec1,... - t,S_max,mu(1),mu(2),sigma(1),sigma(2)); -end - -% define data format -formatSpec2 = '%.4f %.4f %+0.4E %0.4E %1.4f \n'; - -% display global extreme values on screen -fprintf('\n\nGlobal maximum (f Omega PerfFunc power K01):\n\n'); - -fprintf(formatSpec2,... - p1_max_overall,p2_max_overall,... - S_max_overall,power_max_overall,K01_max_overall); - -% display global extreme values in a data file -fprintf(fileID,'\n\nGlobal maximum (f Omega PerfFunc power K01):\n\n'); -fprintf(fileID,formatSpec2,... - p1_max_overall,p2_max_overall,... - S_max_overall,power_max_overall,K01_max_overall); - -fprintf('\n'); - -% close data file -fclose(fileID); - -disp(' '); -time_elapsed = toc -%--------------------------------------------------------------- - - -% save simulation results -% ----------------------------------------------------------- -tic - -disp(' ') -disp(' --- saving simulation results --- '); -disp(' '); -disp(' ... '); -disp(' '); - -save([num2str(case_name),'.mat']); - -toc -% ----------------------------------------------------------- - - - - -% animate cross-entropy iteration -% ........................................................... -%global max: 0.0994 0.7771 +3.4378E-01 3.4378E-01 0.0393 -%global min: 0.0802 0.7771 -8.9856E+00 6.8530E-03 0.9992 -vtitle = 'cross-entropy method'; -%xlab = 'excitation amplitude'; -%ylab = 'excitation frequency'; -xlab = ' mechanical coupling'; -ylab = ' electrical coupling'; -xmin = p1_min; -xmax = p1_max; -ymin = p2_min; -ymax = p2_max; -xref = 0.0994; -yref = 0.7771; -vname = [num2str(case_name),'__animation']; - -plot_ce_animation(p1,p2,(1:t),xref,yref,... - vtitle,xlab,ylab,xmin,xmax,ymin,ymax); -%mov = plot_ce_animation_video(p1,p2,(1:t),... -% xref,yref,vtitle,vname,... -% xlab,ylab,xmin,xmax,ymin,ymax); -% ........................................................... - - - diff --git a/STONEHENGE-1.0/NonconvexOptimization/main_piezomagbeam_opt_ds.m b/STONEHENGE-1.0/NonconvexOptimization/main_piezomagbeam_opt_ds.m deleted file mode 100644 index 760c37a..0000000 --- a/STONEHENGE-1.0/NonconvexOptimization/main_piezomagbeam_opt_ds.m +++ /dev/null @@ -1,421 +0,0 @@ - -% ----------------------------------------------------------------- -% main_piezomagbeam_opt_ds.m -% -% This is the main file for a program which solve an optimization -% problem that seaks to maximize the output power of a -% piezo-magneto-elastic beam, which evolves according to the -% follwing system of ordinary differential equations -% -% d2x/dt2 + 2*ksi*dx/dt - 0.5*x*(1-x^2) - chi*v = f*cos(Omega*t) -% -% dv/dt + lambda*v + kappa*dx/dt = 0 -% -% + -% -% initial conditions, -% -% where -% -% x(t) - dimensionless displacement of the beam tip -% v(t) - dimensionless voltage across the load resistance -% t - dimensionless time -% ksi - mechanical damping ratio -% chi - dimensionless piezoeletric coupling term (mechanical) -% f - dimensionless excitation amplitude -% Omega - dimensionless excitation frequency -% lambda - dimensionless time constant reciprocal -% kappa - dimensionless piezoeletric coupling term (eletrical) -% -% Reference: -% -% A. Erturk, J. Hoffmann, and D. J. Inman -% A piezomagnetoelastic structure for broadband vibration -% energy harvesting -% Applied Physics Letters -% vol. 94 pp. 254102, 2009 -% ----------------------------------------------------------------- -% programmer: Americo Barbosa da Cunha Junior -% americo.cunhajr@gmail.com -% -% last update: Dec 8, 2016 -% ----------------------------------------------------------------- - - -clc -clear -close all - - -% program header -% ----------------------------------------------------------- -disp(' ') -disp(' Piezo-Magneto-Elastic Beam Dynamics ') -disp(' (optimization problem) ') -disp(' ') -disp(' by ') -disp(' Americo Barbosa da Cunha Junior ') -disp(' americo.cunhajr@gmail.com ') -disp(' ') -% ----------------------------------------------------------- - - -% define physical parameters -% ----------------------------------------------------------- -disp(' '); -disp(' --- defining physical parameters --- '); -disp(' '); -disp(' ... '); -disp(' '); - -ksi = 0.01; % mechanical damping ratio -chi = 0.05; % dimensionless piezoeletric coupling term (mechanical) -lambda = 0.05; % dimensionless time constant reciprocal -kappa = 0.5; % dimensionless piezoeletric coupling term (eletrical) -f = 0.0; % dimensionless excitation amplitude -Omega = 0.0; % dimensionless excitation frequency - -x0 = 1.0; % dimensionless initial displacement -xdot0 = 0.0; % dimensionless initial velocity -v0 = 0.0; % dimensionless initial voltage -%--------------------------------------------------------------- - - -% 0-1 test for chaos parameters -% ----------------------------------------------------------- - -% tolerance for 0-1 test -tol01 = 0.1; - -% parameter c for 0-1 test -cmin = 0.0; -cmax = 2*pi; - -% number of test repeats -Nc = 100; - -% oversampling flag -OSflag = 0; -% ----------------------------------------------------------- - - -% define ODE solver parameters -% ----------------------------------------------------------- - -% ODE solver optional parameters -%opt = odeset('RelTol',1.0e-9,'AbsTol',1.0e-6); - -% inicial condition vector -IC = [x0; xdot0; v0]; - -% initial time of analysis -t0 = 0.0; - -% final time of analysis -t1 = 2500.0; - -% time interval of analysis -tspan = [t0 t1]; - -% number of dimensionless time steps -%Ndt = length(tspan); - -% number of steps for steady state -%Nss = round(0.5*Ndt); - -% number of jumps -%Njump = round((Ndt-Nss)/Ndt01); -% ----------------------------------------------------------- - - - -% solve optimization problem via direct search -%--------------------------------------------------------------- -tic - -disp(' '); -disp(' --- optimization problem via direct search --- '); -disp(' '); -disp(' ... '); -disp(' '); - - -% RNG seed -rng_stream = RandStream('mt19937ar','Seed',30081984); -RandStream.setDefaultStream(rng_stream); % Matlab 2009 -%RandStream.setGlobalStream(rng_stream); % Matlab 2013 - -% number of points -Np1 = 256; -Np2 = 256; - -% penalty parameter vector (K01 p1max p1min p2max p2min) -Hpenalty = 10.0; - -% control parameter 1 -%p1_min = 0.08; %camplitude f -%p1_max = 0.10; % amplitude f -p1_min = 0.04; % coupling chi -p1_max = 0.06; % coupling chi - -% control parameter 2 -%p2_min = 0.75; % frequency w -%p2_max = 0.85; % frequency w -p2_min = 0.4; % coupling kappa -p2_max = 0.6; % coupling kappa - -% vector of p1 samples (excitation amplitude) -p1 = linspace(p1_min,p1_max,Np1); - -% vector of p2 samples (excitation frequency) -p2 = linspace(p2_min,p2_max,Np2); - -% preallocate memory for PerfFunc -% (remember: x = columns and y = lines) -S = zeros(Np2,Np1); - -% preallocate memory for power function -% (remember: x = columns and y = lines) -power = zeros(Np2,Np1); - -% preallocate memory for 0-1 classifier -% (remember: x = columns and y = lines) -K01 = zeros(Np2,Np1); - -% initialize PerfFunc maximum/mininum value -S_max_overall = -Inf; -S_min_overall = +Inf; - -% initialize power maximum/minimum value -power_max_overall = -Inf; -power_min_overall = +Inf; - -% initialize K01 maximum/minimum value -K01_max_overall = -Inf; -K01_min_overall = +Inf; - -% initialize power maximum/minimum points -p1_max_overall = 0.0; -p2_max_overall = 0.0; -p1_min_overall = 0.0; -p2_min_overall = 0.0; - - -% simulation information -%case_name = ['piezomagbeam_opt_ds_',... -% 'penal_',num2str(Hpenalty),'_',... -% 'Np1_',num2str(Np1),'_',... -% 'Np2_',num2str(Np2)]; - -case_name = ['piezomagbeam_optdesign_ds_',... - 'penal_',num2str(Hpenalty),'_',... - 'Np1_',num2str(Np1),'_',... - 'Np2_',num2str(Np2)]; - - -% define data file name -file_name = [case_name,'.dat']; - -% open data file -fileID = fopen(file_name,'w'); - -% define data format -formatSpec = '%.4f %.4f %+0.4E %0.4E %1.4f \n'; - -% display extreme values on screen -%fprintf('\n\nLocal values: f Omega PerfFunc power K01\n\n'); - -for np1 = 1:Np1 - - % print loop indicador - disp(['Np1 = ',num2str(np1)]); - - for np2 = 1:Np2 - - % update parameters - %f = p1(np1); - %Omega = p2(np2); - chi = p1(np1); - kappa = p2(np2); - - % define physical paramters vector - phys_param = [ksi chi f Omega lambda kappa]; - - % penalized performance function S(x) - [S(np2,np1),power(np2,np1),K01(np2,np1)] = ... - piezomagbeam_opt_PerfFunc(phys_param,tspan,IC,... - cmin,cmax,Nc,tol01,OSflag,Hpenalty,... - p1_min,p1_max,p2_min,p2_max); - - % update global maximum - if S(np2,np1) > S_max_overall - - S_max_overall = S(np2,np1); - power_max_overall = power(np2,np1); - K01_max_overall = K01(np2,np1); - p1_max_overall = p1(np1); - p2_max_overall = p2(np2); - end - - % update global minimum - if S(np2,np1) < S_min_overall - - S_min_overall = S(np2,np1); - power_min_overall = power(np2,np1); - K01_min_overall = K01(np2,np1); - p1_min_overall = p1(np1); - p2_min_overall = p2(np2); - end - - % print local values in data file - fprintf(fileID,formatSpec,... - p1(np1),p2(np2),S(np2,np1),power(np2,np1),K01(np2,np1)); - - end -end - -% print global maximum -fprintf(fileID,formatSpec,... - p1_max_overall,p2_max_overall,... - S_max_overall,power_max_overall,K01_max_overall); - -% print global minimum -fprintf(fileID,formatSpec,... - p1_min_overall,p2_min_overall,... - S_min_overall,power_min_overall,K01_min_overall); - -% close data file -fclose(fileID); - - -% display global extreme values on screen -fprintf('\n\nGlobal maximum (f Omega PerfFunc power K01):\n\n'); - -fprintf(formatSpec,... - p1_max_overall,p2_max_overall,... - S_max_overall,power_max_overall,K01_max_overall); - -fprintf('\n\nGlobal minimum (f Omega PerfFunc power K01):\n\n'); - -fprintf(formatSpec,... - p1_min_overall,p2_min_overall,... - S_min_overall,power_min_overall,K01_min_overall); - -fprintf('\n'); - -disp(' '); -time_elapsed = toc -%--------------------------------------------------------------- - - - -% save simulation results -% ----------------------------------------------------------- -tic - -disp(' ') -disp(' --- saving simulation results --- '); -disp(' '); -disp(' ... '); -disp(' '); - -save([num2str(case_name),'.mat']); - -toc -% ----------------------------------------------------------- - - -% plot mean power countourmap -% ........................................................... -gtitle = ' performance function contour map'; -%xlab = ' excitation amplitude'; -%ylab = ' excitation frequency'; -xlab = ' mechanical coupling'; -ylab = ' electrical coupling'; -xmin = p1_min; -xmax = p1_max; -ymin = p2_min; -ymax = p2_max; -gname = [case_name,'__perf_func']; -flag = 'eps'; -fig1 = graph_contour_pnt(p1,p2,S',... - p1_max_overall,p2_max_overall,... - gtitle,xlab,ylab,xmin,xmax,ymin,ymax,gname,flag); -%close(fig1); -% ........................................................... - - -% plot mean power countourmap -% ........................................................... -gtitle = ' mean power contour map'; -%xlab = ' excitation amplitude'; -%ylab = ' excitation frequency'; -xlab = ' mechanical coupling'; -ylab = ' electrical coupling'; -xmin = p1_min; -xmax = p1_max; -ymin = p2_min; -ymax = p2_max; -gname = [case_name,'__mean_power']; -flag = 'eps'; -%fig2 = graph_contour_pnt(p1,p2,power',... -% p1_max_overall,p2_max_overall,... -% gtitle,xlab,ylab,xmin,xmax,ymin,ymax,gname,flag); -fig2 = graph_contour_pnt(p1,p2,power',... - p1_max_overall,p2_max_overall,... - gtitle,xlab,ylab,xmin,xmax,ymin,ymax,gname,flag); - -%close(fig2); -% ........................................................... - - -% plot 0-1 identification parameter -% ........................................................... -%xlab = 'excitation amplitude'; -%ylab = 'excitation frequency'; -xlab = ' mechanical coupling'; -ylab = ' electrical coupling'; -label0 = 'regular'; -label1 = 'chaos'; -gtitle = '0-1 test for chaos'; -xmin = p1_min; -xmax = p1_max; -ymin = p2_min; -ymax = p2_max; -zmin = 0; -zmax = 1; -gname = [case_name,'__01test']; -flag = 'eps'; -fig3 = graph_binarymap(p1,p2,K01,... - p1_max_overall,p2_max_overall,... - gtitle,xlab,ylab,label0,label1,... - xmin,xmax,ymin,ymax,zmin,zmax,gname,flag); -% fig3 = graph_binarymap(p1,p2,K01,... -% [],[],... -% gtitle,xlab,ylab,label0,label1,... -% xmin,xmax,ymin,ymax,zmin,zmax,gname,flag); -%close(fig3); -% ........................................................... - - -% plot 0-1 identification parameter ZOOM -% ........................................................... -%xlab = 'excitation amplitude'; -%ylab = 'excitation frequency'; -xlab = ' mechanical coupling'; -ylab = ' electrical coupling'; -label0 = 'regular'; -label1 = 'chaos'; -gtitle = '0-1 test for chaos'; -xmin = 0.099; -xmax = 0.1; -ymin = 0.77; -ymax = 0.78; -gname = [case_name,'__01test_zoom']; -flag = 'eps'; -fig4 = graph_binarymap(p1,p2,K01,... - p1_max_overall,p2_max_overall,... - gtitle,xlab,ylab,label0,label1,... - xmin,xmax,ymin,ymax,zmin,zmax,gname,flag); -%close(fig4); -% ........................................................... diff --git a/STONEHENGE-1.0/NonconvexOptimization/piezomagbeam_opt_PerfFunc.m b/STONEHENGE-1.0/NonconvexOptimization/piezomagbeam_opt_PerfFunc.m deleted file mode 100644 index 7242df7..0000000 --- a/STONEHENGE-1.0/NonconvexOptimization/piezomagbeam_opt_PerfFunc.m +++ /dev/null @@ -1,63 +0,0 @@ - -% ----------------------------------------------------------------- -% piezomagbeam_opt_PerfFunc.m -% -% This function computes the performance function for -% cross entropy method. -% ----------------------------------------------------------------- -% programmer: Americo Barbosa da Cunha Junior -% americo.cunhajr@gmail.com -% -% last update: March 8, 2017 -% ----------------------------------------------------------------- - -% ----------------------------------------------------------------- -function [S,power,K01] = piezomagbeam_opt_PerfFunc(phys_param,... - tspan,IC,... - cmin,cmax,Nc,tol01,... - OSflag,Hpenalty,... - p1_min,p1_max,... - p2_min,p2_max) - -% physical parameters -%ksi = phys_param(1); -%chi = phys_param(2); -f = phys_param(3); -Omega = phys_param(4); -%lambda = phys_param(5); -%kappa = phys_param(6); -%x0 = phys_param(7); -%xdot0 = phys_param(8); -%v0 = phys_param(9); - -% integrate the dynamical system with RK-45 -[time,y] = ode45(@(t,y)piezomagbeam(t,y,phys_param),tspan,IC); - -% number of time steps -Ndt = length(time); - -% number of steps for steady state -Nss = round(0.5*Ndt); - -% number of jumps -Njump = round((Ndt-Nss)/(0.05*Ndt)); - -% temporal interval of analysis -T = time(end)-time(Nss); - -% compute the output power -power = (1/T)*trapz(time(Nss:end),y(Nss:end,3).^2); - -% apply 0-1 test for chaos in voltage time series -K01 = test01chaos(y(Nss:Njump:end,3),cmin,cmax,Nc,OSflag); - -% penalized performance function -S = power - ... - Hpenalty*max(( K01 - tol01),0) + ... - Hpenalty*max(( f - p1_max),0) + ... - Hpenalty*max(( -f + p1_min),0) + ... - Hpenalty*max(( Omega - p2_max),0) + ... - Hpenalty*max((-Omega + p2_min),0); - -end -% ----------------------------------------------------------------- diff --git a/STONEHENGE-1.0/NonconvexOptimization/plot_ce_animation_video.m b/STONEHENGE-1.0/NonconvexOptimization/plot_ce_animation_video.m deleted file mode 100644 index 8b87540..0000000 --- a/STONEHENGE-1.0/NonconvexOptimization/plot_ce_animation_video.m +++ /dev/null @@ -1,121 +0,0 @@ - -% ----------------------------------------------------------------- -% plot_ce_animation_video.m -% -% This functions plots an animation of cross-entropy method -% for optimization -% -% input: -% x1 - x data vector -% y1 - y data vector -% itervec - iteration index vector -% xref - reference value for x -% yref - reference value for y -% vtitle - video title -% vname - video file name -% xlab - x axis label -% ylab - y axis label -% xmin - x axis minimum value -% xmax - x axis maximum value -% ymin - y axis minimum value -% ymax - y axis maximum value -% ----------------------------------------------------------------- -% programmer: Americo Barbosa da Cunha Junior -% americo.cunhajr@gmail.com -% -% last update: Feb 10, 2016 -% ----------------------------------------------------------------- - -% ----------------------------------------------------------------- -function mov = plot_ce_animation_video(x1,y1,itervec,... - xref,yref,vtitle,vname,... - xlab,ylab,xmin,xmax,ymin,ymax) - - % check number of arguments - if nargin < 13 - error('Too few inputs.') - elseif nargin > 13 - error('Too many inputs.') - end - - % check arguments - if size(x1) ~= size(y1) - error('x1 and y1 must be same dimensions') - end - - % coordinates of legend - xiter = xmin + 0.1*(xmax-xmin); - yiter = ymin + 0.9*(ymax-ymin); - - % number of iterations - Niter = length(itervec); - - % open video object - writerObj = VideoWriter(vname); - writerObj.FrameRate = 1; - open(writerObj); - - % preallocate memory for movie frames - nFrames = length(1:Niter); - mov(1:nFrames) = struct('cdata',[], 'colormap',[]); - - % loop to construct the video - for n=1:Niter - - % initialize video frame - fig = figure(1000); - - % define frame properties - set(gcf,'color','white'); - set(gca,'XColor',[.3 .3 .3],'YColor',[.3 .3 .3],'zColor',[.3 .3 .3]); - set(gca,'FontName','Helvetica'); - set(gca,'FontSize',16); - set(gca,'Box','on'); - - % draw trailer base right arm - fh0 = plot(x1(:,n),y1(:,n),'xb'); - - hold on - - fh1 = plot(xref,yref,'xr'); - - set(fh0,'LineWidth',2.0); - set(fh0,'MarkerSize',10.0); - set(fh0,'MarkerFaceColor','w'); - set(fh0,'MarkerEdgeColor','b'); - set(fh1,'LineWidth',5.0); - set(fh1,'MarkerSize',20.0); - set(fh1,'MarkerFaceColor','w'); - set(fh1,'MarkerEdgeColor','r'); - - % define labels - xlabel(xlab,'FontSize',22,'FontName','Helvetica'); - ylabel(ylab,'FontSize',22,'FontName','Helvetica'); - - % define video frame limits - xlim([xmin xmax]); - ylim([ymin ymax]); - - % display time counter - text(xiter,yiter,['iteration = ',num2str(itervec(n))],... - 'Color','k',... - 'FontName','Helvetica',... - 'FontSize',12); - - hold off - - % define video title - title(vtitle,'FontSize',20,'FontName','Helvetica'); - - % save movie frame - mov(n) = getframe(gcf); - writeVideo(writerObj,mov(n)); - - pause(1.5) - end - - % close movie object - close(writerObj); -return -% ----------------------------------------------------------------- -