基于MeSC与交感神经作用关系的压力水平与白发模拟系统和压力规划系统（Matlab）
作者：互联网
基于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