标签：nh ... whr fix bi hair MeSC 模拟系统 规划系统
基于MeSC与交感神经作用关系的压力水平与白发模拟系统和压力规划系统
一个本人的Matlab项目，可用于根据压力水平模拟白发水平，并根据工作情况给出白发量最少的合理的压力规划。
% 细胞仿真
clc;
clear;
% Raw data
sl0 = [0 1 2 3 4]; % stress_level
MeSC = [7.5 10 5 2.5 1]; %No. of MeScs within per HF
whr1 = 0.4; % max white hair rate 已知MeSC为1时白发率为40%
whr0 = 0; % whr when no stress 设无压力时白发率为0
Fr = 0.33; % 女性的白发率
Mr = 0.29; % 男性的白发率
% spline插值获得白发率
x = [MeSC(end) MeSC(1)];
y = [whr1 whr0];
whr = spline(x,y,MeSC); % 利用已知的白发率与MeSC的关系插值
% Data visualization
figure()
y = [sl0;MeSC];
X = categorical({'Group1','Group2','Group3','Group4','Group5'});
b = bar(X,y);hold on
xtips1 = b(1).XEndPoints;
ytips1 = b(1).YEndPoints;
labels1 = string(b(1).YData);
text(xtips1,ytips1,labels1,'HorizontalAlignment','center',...
'VerticalAlignment','bottom')
xtips2 = b(2).XEndPoints;
ytips2 = b(2).YEndPoints;
labels2 = string(b(2).YData);
text(xtips2,ytips2,labels2,'HorizontalAlignment','center',...
'VerticalAlignment','bottom')
ylim([0,12]);
grid on;
title('Experiment data');
xlabel('test group');
colororder({'k','k'})
yyaxis left
ylabel('Stress level & MeSC level');
yyaxis right
ylabel('White hair rate');
plot(X,whr,'linewidth',2,'color','b');hold on; % 利用压力水平与MeSC的对应关系，得到压力水平和白发率关系
scatter(X,whr);hold on;
ylim([-0.2,0.5]);
legend('stress level','MeSC level','White hair rate');
% Parameter setting
upr = 0.0001; % 头发更新率 ???
nh = 100000; % 头发数量，一般10万
nh_l = sqrt(nh); % 生成图片的长度
Initial = 20;
%{
% Stress data input & simulation
% sli = input('每周的压力水平？(0-4)\n');
sli = input('每周的睡眠分数？(0-4)\n');
sli = ((100-sli)./100).*4;
sli = sli';
n = size(sli);
n0 = 1:1:n; % 已知点
n1 = 1:n/(7*n(1)):n; % 插值点
sl = spline(n0,sli,n1); % 对压力水平插值
nsl = size(sl);
% Stress visualization
figure();
grid on;
yyaxis left
p1 = plot(1:nsl(2),sl,'linewidth',1.25,'color','b'); % 每日压力情况曲线
hold on;
stem(1:7:nsl(2),sli,'linewidth',1.25,'LineStyle','-.',...
'MarkerFaceColor','red',...
'MarkerEdgeColor','green');hold on;
title('Personal data');
xlabel('days');
ylabel('stress level');
% White hair rate simulation & visualization
whr2 = spline(sl0,whr,sl); % 根据输入的压力情况得到实际的每日白发率曲线
nwhr2 = size(whr2);
yyaxis right
p2 = plot(1:nsl(2),whr2,'linewidth',1.5,'color','r');hold on;
ylabel('white hair rate');
colororder({'b','r'})
legend([p1,p2],{'stress level','White hair rate','MeSC level'});
MeSC2 = spline(whr,MeSC,whr2); % 插值估计每日MeSC数量
% p3 = plot(1:nsl(2),MeSC2,'linewidth',1.5,'color','k');hold on;
%{
% 椭球头部模型
figure()
[x, y, z] = ellipsoid(0,0,0,3,4,4.5,fix(nh_l));
surf(x, y, -z)
c = copper;
c = flipud(c);
colormap(c);
axis equal
%}
% Hair simulation & visualization
hair = zeros(fix(nh_l),fix(nh_l))+Initial; % fix取整,设置头发初始灰度
hair = uint8(hair);
wi = 1; % white i,变白了的头发的索引
bi = 1; % black i,变黑了的头发的索引
%F = zeros(1,nsl(2)*fix(upr*nh));
F = moviein(nsl(2)*fix(upr*nh));
j = 1;
wx = randi(fix(nh_l),1,fix(nh_l));
wy = randi(fix(nh_l),1,fix(nh_l));
for t = 1:1:nsl(2) % 对每一天
x = fix((nh_l-1)*rand(upr*nh,1)+1); % 求该天被更新的头发坐标集合，rand生成数组
y = fix((nh_l-1)*rand(upr*nh,1)+1);
for i = 1:fix(upr*nh)
i_whr = whr2(t);
if i_whr > 0
hair(x(i),y(i)) = hair(x(i),y(i))+255*i_whr;
wx(wi) = x(i); % 记录下变白了的头发
wy(wi) = y(i);
wi = wi+1;
else
if (hair(wx(bi),wy(bi))+255*i_whr)>=Initial
hair(wx(bi),wy(bi)) = ...
hair(wx(bi),wy(bi))+255*i_whr; % 此前变白的头发将优先恢复
bi = bi+1;
end
end
hair2 = mat2gray(hair);
fig = figure();
fig.Visible = 'off';
subplot(121)
imshow(hair2);
subplot(122)
Mf = imread('MeSC.png');
imshow(Mf);
hold on;
N = MeSC2(t);
ym = 829*(1-0.45).*rand(round(N),1);
xm = ones(1,round(N)).*0.7*318;
xmr= -20+(20+20)*rand(1,round(N));
xm = xm+xmr;
scatter(xm,ym,60,'g','filled'); % 模拟MeSC的数量与分布
F(j) = getframe(fig);
j = j+1;
end
end
figure()
%{
set(gcf,'unit','normalized',... % 本来打算设置figure大小
'position',[0.1,0.1,0.56,0.42])
xlim([0,0.5]);
ylim([0,0.8]);
%}
movie(F,5);
%}
% Optimization
% 设对每一时刻，一个迫近（一周内）ddl的压力值为0.7，第二周的ddl压力值为0.3
% 对未来两周
% ddl input & stress level count
ddl1 = input('第一周ddl数目：');
ddl2 = input('第二周ddl数目：');
ddl3 = input('第三周ddl数目：');
sl1 = 0.7*ddl1+0.3*ddl2;
sl2 = 0.7*ddl2+0.3*ddl3;
% 多项式拟合求出sl-whr函数表达式
p=polyfit(sl0,whr,4);
whr1 = polyval(p,sl0);
figure()
plot(sl0,whr1,'linewidth',5,'color','b');hold on
plot(sl0,whr,'linewidth',2.5,'color','g');grid on
title('Polynomial fitting');
legend('Polynomial output','truth');
xlabel('stress level');
ylabel('white hair rate');
fun = @(x)-(7*p(1)+p(2)*(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)+x(9)+x(10)+x(11)+x(12)+x(13)+x(14))...
+p(2)*(x(1)^2+x(2)^2+x(3)^2+x(4)^2+x(5)^2+x(6)^2+x(7)^2+x(8)^2+x(9)^2+x(10)^2+x(11)^2+x(12)^2+x(13)^2+x(14)^2)...
+p(3)*(x(1)^3+x(2)^3+x(3)^3+x(4)^3+x(5)^3+x(6)^3+x(7)^3+x(8)^3+x(9)^3+x(10)^3+x(11)^3+x(12)^3+x(13)^3+x(14)^3)...
+p(4)*(x(1)^4+x(2)^4+x(3)^4+x(4)^4+x(5)^4+x(6)^4+x(7)^4+x(8)^4+x(9)^4+x(10)^4+x(11)^4+x(12)^4+x(13)^4+x(14)^4));
x0 = [0,0,0,0,0,0,0,0,0,0,0,0,0,0];
A = [-1,-1,-1,-1,-1,-1,-1,0,0,0,0,0,0,0;... % 保证第一周ddl完成
0,0,0,0,0,0,0,-1,-1,-1,-1,-1,-1,-1;... % 保证第二周ddl完成
1,1,0,0,0,0,0,0,0,0,0,0,0,0;... % 别太累了，连续两天压力值不超过6
0,1,1,0,0,0,0,0,0,0,0,0,0,0;...
0,0,1,1,0,0,0,0,0,0,0,0,0,0;...
0,0,0,1,1,0,0,0,0,0,0,0,0,0;...
0,0,0,0,1,1,0,0,0,0,0,0,0,0;...
0,0,0,0,0,1,1,0,0,0,0,0,0,0;...
0,0,0,0,0,0,1,1,0,0,0,0,0,0;...
0,0,0,0,0,0,0,1,1,0,0,0,0,0;...
0,0,0,0,0,0,0,0,1,1,0,0,0,0;...
0,0,0,0,0,0,0,0,0,1,1,0,0,0;...
0,0,0,0,0,0,0,0,0,0,1,1,0,0;...
0,0,0,0,0,0,0,0,0,0,0,1,1,0;...
0,0,0,0,0,0,0,0,0,0,0,0,1,1;]; % 不等式约束的系数矩阵，>号，要取相反数
b = [-sl1*7;-sl2*7;6;6;6;6;6;6;6;6;6;6;6;6;6]; % 不等式约束的常向量，>号，要取相反数
Aeq = []; % 等式约束的系数矩阵
beq = []; % 等式约束的常向量
lb = [0;0;0;0;0;0;0;0;0;0;0;0;0;0;]; % 自变量的最小值
ub = [4;4;4;4;4;4;4;4;4;4;4;4;4;4;]; % 自变量的最大值
[x,fval] = fmincon(fun,x0,A,b,Aeq,beq,lb,ub);
figure();
subplot(2,2,[1,2])
yyaxis left
plot(1:14,x,'linewidth',1.5,'color','b');hold on; % 最合理的压力规划
grid on;
ylabel('stress level');
xlabel('days');
yyaxis right
whr_c = zeros(1,14);
for i = 1:14
for k = 1:i
whr_c(i) = whr_c(i)+polyval(p,x(k));
end
end
plot(1:14,whr_c*upr,'linewidth',1.5);hold on; % 累计白发率
ylabel('cumulative white hair rate');
title('Optimization result');
% hair simulation VS
hair = zeros(fix(nh_l),fix(nh_l))+Initial; % fix取整,设置头发初始灰度
hair = uint8(hair);
wi = 1; % white i,变白了的头发的索引
bi = 1; % black i,变黑了的头发的索引
wx = randi(fix(nh_l),1,fix(nh_l));
wy = randi(fix(nh_l),1,fix(nh_l));
for t = 1:1:14 % 对每一天
xx = fix((nh_l-1)*rand(upr*nh,1)+1); % 求该天被更新的头发坐标集合，rand生成数组
yy = fix((nh_l-1)*rand(upr*nh,1)+1);
for i = 1:fix(upr*nh)
i_whr = polyval(p,x(t));
if i_whr > 0
hair(xx(i),yy(i)) = hair(xx(i),yy(i))+255*i_whr;
wx(wi) = xx(i); % 记录下变白了的头发
wy(wi) = yy(i);
wi = wi+1;
else
hair(wx(bi),wy(bi)) = ...
hair(wx(bi),wy(bi))+255*i_whr; % 此前变白的头发将优先恢复
bi = bi+1;
end
end
end
subplot(2,2,3)
imshow(hair);
xlabel('optimal result');
hair = zeros(fix(nh_l),fix(nh_l))+Initial; % fix取整,设置头发初始灰度
hair = uint8(hair);
wi = 1; % white i,变白了的头发的索引
bi = 1; % black i,变黑了的头发的索引
wx = randi(fix(nh_l),1,fix(nh_l));
wy = randi(fix(nh_l),1,fix(nh_l));
for t = 1:1:14 % 对每一天
xx = fix((nh_l-1)*rand(upr*nh,1)+1); % 求该天被更新的头发坐标集合，rand生成数组
yy = fix((nh_l-1)*rand(upr*nh,1)+1);
for i = 1:fix(upr*nh)
i_whr = (sl1+sl2)/2;
if i_whr > 0
hair(xx(i),yy(i)) = hair(xx(i),yy(i))+255*i_whr;
wx(wi) = xx(i); % 记录下变白了的头发
wy(wi) = yy(i);
wi = wi+1;
else
hair(wx(bi),wy(bi)) = ...
hair(wx(bi),wy(bi))+255*i_whr; % 此前变白的头发将优先恢复
bi = bi+1;
end
end
end
subplot(2,2,4)
imshow(hair);
xlabel('Non-optimal result result');
标签：nh,...,whr,fix,bi,hair,MeSC,模拟系统,规划系统 来源： https://blog.csdn.net/weixin_46813313/article/details/118667823
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。