思维创新策略（Thinking Innovation Strategy, TIS） 源于对人类思维方式的深刻洞察。面对复杂而多变的问题，人类往往并非一味追求唯一的“完美答案”，而是善于借鉴成功经验，汲取他人的有效策略，通过不断筛选与调整，在全局中寻找更优的可能。TIS 正是将这一“个体思维—集体进化”的过程抽象为算法机制：不直接锁定最优解，而是以表现出色的个体为核心，通过“适者生存”推动解集向更优方向演化。与传统依赖探索与开发平衡的元启发式方法不同，TIS 以思维的创新性驱动优化过程，让算法在面对高复杂度、多约束的工程问题时，依旧能够展现出稳健的全局搜索能力与持续的性能提升潜力。
该成果于2025年7月发表于计算机领域二区sci期刊Applied Soft Computing，具有广阔的应用前景。
TIS 是一种基于人类思维现象的启发式策略，尽管已有许多改进型通用策略，但 TIS 在本质上与这些策略存在显著差异，并能够帮助算法找到更优的解。人类的认知过程与思维模式十分复杂，涵盖感知、知觉、注意、记忆、思考等多个方面。此外，人类思维模式中还存在多种知识深度（Depth of Knowledge, DOK）层次，从回答仅需有限技能与知识的简单问题，到处理复杂且新颖的问题不等。从神经科学的角度来看，人类思维的本质可通过神经元之间的电信号传递及其连接模式来解释，这些连接共同构成了大脑中的认知与思维框架。
1、算法原理
（1）创新思维
创新思维（又称“创造性思维”）指的是发明或发现新方法、新表达的过程。它涉及将对创新意识的渴望从感性探索提升至理性探究，推动认知由感性阶段向理性阶段转化。这一高级思维形式是人类探索现象本质、获取新知识、发展新能力的有效途径。当面对多样且复杂的问题时，个体常依赖其知识深度（Depth of Knowledge, DOK）来寻找解决方案。DOK 的作用相当于一个支持批判性思维的综合性心理资源库，其主要包含两个层次：DOK1：个人成长过程中积累的知识与经验。DOK2：关于自然现象的持续获取与信息储存。
这种结构化表征突出了 DOK 如何凭借其强大的认知能力，使个体能够有效应对传统的简单问题以及复杂的新颖挑战。DOK 的计算方式如下所示：
其中，DOK1 表示随时间积累的知识与经验；DOK2 表示持续获取与存储的信息； 表示常数 ；FEs 表示当前评估次数，MaxFEs 表示最大评估次数。想象力（IM）指通过新颖方式协调感知材料以创造新形象的能力。个体想象力对于创新思维至关重要，其发展源于对既有知识与经验的批判性评价，以及对研究对象的全面观察与理解。在个体认知的基础上，它通过直觉、表象或多维度结合进行整合，不受既有结论、观点、框架或理论的限制，提出新的构想。IM 的计算公式如下式所示：
其中，rand 表示 0 到 1 之间的随机数。创新性思维源于想象力、遗传禀赋、天赋、智商及思维方式之间的相互作用，并逐渐发展为一种强调新颖性、发散性与适切性的原创思维形式。批判性思维的培养能够促使个体质疑现有信息，从而激发新的科学认知想象与灵感。当个体旨在产出具有新颖性与价值的成果时，常会将其深厚的知识积累与丰富的想象力相结合或加以转化，以实现创新。创新性思维的计算公式如下式所示：
其中， 表示更新后种群的位置， 表示当前个体的位置。
（2）个体思维过程展示
下图展示了一个模拟个体解决问题的思维步骤。在个体思维阶段，大脑首先读取与问题相关的信息并保留信息事件。随后，计算个体成长过程中所吸收的经验知识及所保留自然信息的知识深度（DOK）。第三，基于当前算法状态，实施创新思维，增强基本算法探索阶段的开发能力。与此同时，在开发阶段，通过创新思维扩大个体搜索范围，并在有效空间内开展多方向探索，此举有助于算法避免陷入局部最优。最终，借助创新思维获得更有效的个体解，从而提升算法的优化性能。总之，引入创新思维有效引导个体进行多向思考，最大化问题空间的探索与开发能力，增强了算法的适应性和实用性，使其适用于多种类型的优化问题。
下图展示了创新思维的创造过程。图中， 代表当前个体的位置， 代表通过思维创新策略获得的最新个体解的位置。“ER” 表示基本算法的探索阶段，“ET” 表示开发阶段，“IT” 则表示个体执行位置更新时创新思维策略的实施。
（3）知识深度层级
个体思维源于人类对外部刺激的认知。常用且简单的知识通常储存在大脑的第一层和第二层知识中，代表通过学习获得的经验与技能。而对于未知或较少接触的信息，如新的自然现象或非常规知识，则通常储存在第三层和第四层知识中。这种更深层次的知识在解决陌生问题时，是创新的重要源泉。个体思维知识的层级结构如下图所示。
（4）TIS 的伪代码及流程图实现
为了更全面地展示思维创新策略（Thinking Innovation Strategy，简称TIS）的应用效果，本案例选取了一种基于TIS的混合进化算法——即将TIS引入混合自适应差分进化与协方差矩阵适应进化策略（TIS_LSHADE_SPACMA）作为示范，详细阐述TIS在优化算法中的应用方法及实现流程。
对于种群中个体 ，通过差分进化产生变异向量，其计算公式为：
其中， 为随机选择的个体索引， 表示从当前前 %优良个体中随机选择的指导个体， 为缩放因子。
同时，算法利用协方差矩阵适应进化策略（CMA-ES）对变异向量进行局部搜索，基于多变量正态分布：
其中， 为步长参数， 和 分别为协方差矩阵的特征向量矩阵与对角特征值矩阵， 为均值向量。
本算法的具体操作流程以伪代码形式呈现（见算法1），涵盖了种群初始化、个体位置更新、适应度评估、存档管理、记忆库更新及动态种群规模调整等关键步骤。值得强调的是，TIS模块在每次个体位置更新时均被调用，通过创新思维驱动对个体搜索位置的多维再探索，有效提升算法跳出局部最优的能力。
为便于理解具体实现过程，基于伪代码绘制了该基于TIS的混合进化算法流程图，系统展示了算法的整体运行机制与核心流程。
2、MATLAB核心代码
.rtcContent { padding: 30px; } .lineNode {font-size: 14pt; font-family: "新宋体", Menlo, Monaco, Consolas, "Courier New", monospace; font-style: normal; font-weight: normal; }function [gbestval,gbest,run_funcvals]= TIS_LSHADE_SPACMA(N,FESmax,lb,ub,dim,fobj)L_Rate= 0.80;lu = [lb; ub];p_best_rate = 0.11;arc_rate = 1.4;memory_size = 5;pop_size = 18 * dim;max_pop_size = pop_size;min_pop_size = 4.0;First_calss_percentage=0.5;popold = repmat(lu(1, :), pop_size, 1) + rand(pop_size, dim) .* (repmat(lu(2, :) - lu(1, :), pop_size, 1));pop = popold;gbestval = inf;gbest = zeros(1, dim);fitness=[];run_funcvals=[];count=0;pmax=FESmax/10;person=[];personf=[];lu=[lb;ub]; %problem boundary information% lu=repmat(lu,1,dim);fitness(1,1)=feval(fobj,pop(1,:));FES=1;m=1;person(m,:)=pop(1,:); %store the first memorypersonf(m,1)=fitness(1,1); for i = 1 : pop_size [m,person,personf,popold(i,:),fitness(i,:),FES]=... TIS(pmax,m,person,personf,popold(i,:),fobj,FES,lu,FESmax); if (fitness(i) < gbestval && isreal(pop(i, :)) && sum(isnan(pop(i, :)))==0 ) gbestval = fitness(i); gbest = popold(i, :); endendrun_funcvals = ones(1,FES-length(run_funcvals))*gbestval;memory_sf = 0.5 .* ones(memory_size, 1);memory_cr = 0.5 .* ones(memory_size, 1);memory_pos = 1;archive.NP = arc_rate * pop_size; % the maximum size of the archivearchive.pop = zeros(0, dim); % the solutions stored in te archivearchive.funvalues = zeros(0, 1); % the function value of the archived solutionsmemory_1st_class_percentage = First_calss_percentage.* ones(memory_size, 1); % Class#1 probability for Hybridization%% Initialize CMAES parameterssigma = 0.5; % coordinate wise standard deviation (step size)xmean = rand(dim,1); % objective variables initial pointmu = pop_size/2; % number of parents/points for recombinationweights = log(mu+1/2)-log(1:mu)'; % muXone array for weighted recombinationmu = floor(mu);weights = weights/sum(weights); % normalize recombination weights arraymueff=sum(weights)^2/sum(weights.^2); % variance-effectiveness of sum w_i x_i% Strategy parameter setting: Adaptationcc = (4 + mueff/dim) / (dim+4 + 2*mueff/dim); % time constant for cumulation for Ccs = (mueff+2) / (dim+mueff+5); % t-const for cumulation for sigma controlc1 = 2 / ((dim+1.3)^2+mueff); % learning rate for rank-one update of Ccmu = min(1-c1, 2 * (mueff-2+1/mueff) / ((dim+2)^2+mueff)); % and for rank-mu updatedamps = 1 + 2*max(0, sqrt((mueff-1)/(dim+1))-1) + cs; % damping for sigma usually close to 1% Initialize dynamic (internal) strategy parameters and constantspc = zeros(dim,1);ps = zeros(dim,1); % evolution paths for C and sigmaB = eye(dim,dim); % B defines the coordinate systemD = ones(dim,1); % diagonal D defines the scalingC = B * diag(D.^2) * B'; % covariance matrix CinvsqrtC = B * diag(D.^-1) * B'; % C^-1/2eigeneval = 0; % track update of B and DchiN=dim^0.5*(1-1/(4*dim)+1/(21*dim^2)); % expectation of%% main loopHybridization_flag=1; % Indicator flag if we need to Activate CMAES Hybridizationwhile 1 pop = popold; % the old population becomes the current population [~, sorted_index] = sort(fitness, 'ascend'); mem_rand_index = ceil(memory_size * rand(pop_size, 1)); mu_sf = memory_sf(mem_rand_index); mu_cr = memory_cr(mem_rand_index); mem_rand_ratio = rand(pop_size, 1); %% for generating crossover rate cr = normrnd(mu_cr, 0.1); term_pos = find(mu_cr == -1); cr(term_pos) = 0; cr = min(cr, 1); cr = max(cr, 0); %% for generating scaling factor if(FES <= FESmax/2) sf=0.45+.1*rand(pop_size, 1); pos = find(sf <= 0); while ~ isempty(pos) sf(pos)=0.45+0.1*rand(length(pos), 1); pos = find(sf <= 0); end else sf = mu_sf + 0.1 * tan(pi * (rand(pop_size, 1) - 0.5)); pos = find(sf <= 0); while ~ isempty(pos) sf(pos) = mu_sf(pos) + 0.1 * tan(pi * (rand(length(pos), 1) - 0.5)); pos = find(sf <= 0); end end sf = min(sf, 1); %% for generating Hybridization Class probability Class_Select_Index=(memory_1st_class_percentage(mem_rand_index)>=mem_rand_ratio); if(Hybridization_flag==0) Class_Select_Index=or(Class_Select_Index,~Class_Select_Index);%All will be in class#1 end r0 = 1 : pop_size; popAll = [pop; archive.pop]; [r1, r2] = gnR1R2(pop_size, size(popAll, 1), r0); pNP = max(round(p_best_rate * pop_size), 2); %% choose at least two best solutions randindex = ceil(rand(1, pop_size) .* pNP); %% select from [1, 2, 3, ..., pNP] randindex = max(1, randindex); %% to avoid the problem that rand = 0 and thus ceil(rand) = 0 pbest = pop(sorted_index(randindex), :); %% randomly choose one of the top 100p% solutions vi=[]; temp=[]; if(sum(Class_Select_Index)~=0) vi(Class_Select_Index,:) = pop(Class_Select_Index,:) + sf(Class_Select_Index, ones(1, dim)) .* (pbest(Class_Select_Index,:) - pop(Class_Select_Index,:) + pop(r1(Class_Select_Index), :) - popAll(r2(Class_Select_Index), :)); end if(sum(~Class_Select_Index)~=0) for k=1:sum(~Class_Select_Index) temp(:,k) = xmean + sigma * B * (D .* randn(dim,1)); % m + sig * Normal(0,C) end vi(~Class_Select_Index,:) = temp'; end if(~isreal(vi)) Hybridization_flag=0; continue; end vi = boundConstraint(vi, pop, lu); mask = rand(pop_size, dim) > cr(:, ones(1, dim)); % mask is used to indicate which elements of ui comes from the parent rows = (1 : pop_size)'; cols = floor(rand(pop_size, 1) * dim)+1; % choose one position where the element of ui doesn't come from the parent jrand = sub2ind([pop_size dim], rows, cols); mask(jrand) = false; ui = vi; ui(mask) = pop(mask); children_fitness=[]; for i = 1 : pop_size [m,person,personf,ui(i,:),children_fitness(i,:),FES]=... TIS(pmax,m,person,personf,ui(i,:),fobj,FES,lu,FESmax); if (children_fitness(i) < gbestval && isreal(ui(i, :)) && sum(isnan(ui(i, :)))==0) gbestval = children_fitness(i); gbest = ui(i, :); end if FES >= FESmax run_funcvals = [run_funcvals,ones(1,FES-length(run_funcvals))*gbestval]; return end end
run_funcvals = [run_funcvals,ones(1,FES-length(run_funcvals))*gbestval]; dif = abs(fitness - children_fitness); Child_is_better_index = (fitness > children_fitness); goodCR = cr(Child_is_better_index == 1); goodF = sf(Child_is_better_index == 1); dif_val = dif(Child_is_better_index == 1); dif_val_Class_1 = dif(and(Child_is_better_index,Class_Select_Index) == 1); dif_val_Class_2 = dif(and(Child_is_better_index,~Class_Select_Index) == 1); archive = updateArchive(archive, popold(Child_is_better_index == 1, :), fitness(Child_is_better_index == 1)); [fitness, Child_is_better_index] = min([fitness, children_fitness], [], 2); popold = pop; popold(Child_is_better_index == 2, :) = ui(Child_is_better_index == 2, :); num_success_params = numel(goodCR); if num_success_params > 0 sum_dif = sum(dif_val); dif_val = dif_val / sum_dif; %% for updating the memory of scaling factor memory_sf(memory_pos) = (dif_val' * (goodF .^ 2)) / (dif_val' * goodF); %% for updating the memory of crossover rate if max(goodCR) == 0 || memory_cr(memory_pos) == -1 memory_cr(memory_pos) = -1; else memory_cr(memory_pos) = (dif_val' * (goodCR .^ 2)) / (dif_val' * goodCR); end if (Hybridization_flag==1)% if the Hybridization is activated memory_1st_class_percentage(memory_pos) = memory_1st_class_percentage(memory_pos)*L_Rate+ (1-L_Rate)*(sum(dif_val_Class_1) / (sum(dif_val_Class_1) + sum(dif_val_Class_2))); memory_1st_class_percentage(memory_pos)=min(memory_1st_class_percentage(memory_pos),0.8); memory_1st_class_percentage(memory_pos)=max(memory_1st_class_percentage(memory_pos),0.2); end memory_pos = memory_pos + 1; if memory_pos > memory_size; memory_pos = 1; end end %% for resizing the population size plan_pop_size = round((((min_pop_size - max_pop_size) / FESmax) * FES) + max_pop_size); if pop_size > plan_pop_size reduction_ind_num = pop_size - plan_pop_size; if pop_size - reduction_ind_num < min_pop_size; reduction_ind_num = pop_size - min_pop_size;end pop_size = pop_size - reduction_ind_num; for r = 1 : reduction_ind_num [~, indBest] = sort(fitness, 'ascend'); worst_ind = indBest(end); popold(worst_ind,:) = []; pop(worst_ind,:) = []; fitness(worst_ind,:) = []; Child_is_better_index(worst_ind,:) = []; end archive.NP = round(arc_rate * pop_size); if size(archive.pop, 1) > archive.NP rndpos = randperm(size(archive.pop, 1)); rndpos = rndpos(1 : archive.NP); archive.pop = archive.pop(rndpos, :); end %% update CMA parameters mu = pop_size/2; % number of parents/points for recombination weights = log(mu+1/2)-log(1:mu)'; % muXone array for weighted recombination mu = floor(mu); weights = weights/sum(weights); % normalize recombination weights array mueff=sum(weights)^2/sum(weights.^2); % variance-effectiveness of sum w_i x_i end %% CMAES Adaptation if(Hybridization_flag==1) % Sort by fitness and compute weighted mean into xmean [~, popindex] = sort(fitness); % minimization xold = xmean; xmean = popold(popindex(1:mu),:)' * weights; % recombination, new mean value % Cumulation: Update evolution paths ps = (1-cs) * ps + sqrt(cs*(2-cs)*mueff) * invsqrtC * (xmean-xold) / sigma; hsig = sum(ps.^2)/(1-(1-cs)^(2*FES/pop_size))/dim < 2 + 4/(dim+1); pc = (1-cc) * pc + hsig * sqrt(cc*(2-cc)*mueff) * (xmean-xold) / sigma; % Adapt covariance matrix C artmp = (1/sigma) * (popold(popindex(1:mu),:)' - repmat(xold,1,mu)); % mu difference vectors C = (1-c1-cmu) * C ... % regard old matrix + c1 * (pc * pc' ... % plus rank one update + (1-hsig) * cc*(2-cc) * C) ... % minor correction if hsig==0 + cmu * artmp * diag(weights) * artmp'; % plus rank mu update % Adapt step size sigma sigma = sigma * exp((cs/damps)*(norm(ps)/chiN - 1)); % Update B and D from C if FES - eigeneval > pop_size/(c1+cmu)/dim/10 % to achieve O(problem_size^2) eigeneval = FES; C = triu(C) + triu(C,1)'; % enforce symmetry if(sum(sum(isnan(C)))>0 || sum(sum(~isfinite(C)))>0 || ~isreal(C)) Hybridization_flag=0; continue; end [B,D] = eig(C); % eigen decomposition, B==normalized eigenvectors D = sqrt(diag(D)); % D contains standard deviations now invsqrtC = B * diag(D.^-1) * B'; end endend disp("count = "+count);endfunction archive = updateArchive(archive, pop, funvalue)% Update the archive with input solutions% Step 1: Add new solution to the archive% Step 2: Remove duplicate elements% Step 3: If necessary, randomly remove some solutions to maintain the archive size%% Version: 1.1 Date: 2008/04/02% Written by Jingqiao Zhang (jingqiao@gmail.com)if archive.NP == 0, return; endif size(pop, 1) ~= size(funvalue,1), error('check it'); end% Method 2: Remove duplicate elementspopAll = [archive.pop; pop ];funvalues = [archive.funvalues; funvalue ];[~, IX]= unique(popAll, 'rows');if length(IX) < size(popAll, 1) % There exist some duplicate solutions popAll = popAll(IX, :); funvalues = funvalues(IX, :);endif size(popAll, 1) <= archive.NP % add all new individuals archive.pop = popAll; archive.funvalues = funvalues;else % randomly remove some solutions rndpos = randperm(size(popAll, 1)); % equivelent to "randperm"; rndpos = rndpos(1 : archive.NP); archive.pop = popAll (rndpos, :); archive.funvalues = funvalues(rndpos, :);endend%This function is used for L-SHADE bound checkingfunction vi = boundConstraint (vi, pop, lu)% if the boundary constraint is violated, set the value to be the middle% of the previous value and the bound[NP, ~] = size(pop); % the population size and the problem's dimension%% check the lower boundxl = repmat(lu(1, :), NP, 1);pos = vi < xl;vi(pos) = (pop(pos) + xl(pos)) / 2;%% check the upper boundxu = repmat(lu(2, :), NP, 1);pos = vi > xu;vi(pos) = (pop(pos) + xu(pos)) / 2;endfunction [r1, r2] = gnR1R2(NP1, NP2, r0)% gnA1A2 generate two column vectors r1 and r2 of size NP1 & NP2, respectively% r1's elements are choosen from {1, 2, ..., NP1} & r1(i) ~= r0(i)% r2's elements are choosen from {1, 2, ..., NP2} & r2(i) ~= r1(i) & r2(i) ~= r0(i)%% Call:% [r1 r2 ...] = gnA1A2(NP1) % r0 is set to be (1:NP1)'% [r1 r2 ...] = gnA1A2(NP1, r0) % r0 should be of length NP1%% Version: 2.1 Date: 2008/07/01% Written by Jingqiao Zhang (jingqiao@gmail.com)NP0 = length(r0);r1 = floor(rand(1, NP0) * NP1) + 1;%for i = 1 : inffor i = 1 : 99999999 pos = (r1 == r0); if sum(pos) == 0 break; else % regenerate r1 if it is equal to r0 r1(pos) = floor(rand(1, sum(pos)) * NP1) + 1; end if i > 1000 % this has never happened so far error('Can not genrate r1 in 1000 iterations'); endendr2 = floor(rand(1, NP0) * NP2) + 1;%for i = 1 : inffor i = 1 : 99999999 pos = ((r2 == r1) | (r2 == r0)); if sum(pos)==0 break; else % regenerate r2 if it is equal to r0 or r1 r2(pos) = floor(rand(1, sum(pos)) * NP2) + 1; end if i > 1000 % this has never happened so far error('Can not genrate r2 in 1000 iterations'); endendend
参考文献
[1] Jia H, Zhou X, Zhang J. Thinking Innovation Strategy (TIS): A novel mechanism for metaheuristic algorithm design and evolutionary update[J]. Applied Soft Computing, 2025, 175: 113071.
完整代码获取
后台回复关键词：
TGDM879
获取更多代码：
扫一扫
1118
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
