D-Lift Problems

本文介绍了一个名为LiftProblems的问题,该问题关注的是如何最小化乘客在使用电梯过程中的不满情绪。通过动态规划的方法,文章详细阐述了状态转移方程的设计与实现,旨在寻找最优的电梯停靠策略。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

Lift Problems

题目大意:每个学生乘坐电梯如果不在自己想要停的楼层停下就会产生怒气值,在低于自己楼层停下,在停下的楼层就会产生一个怒气值,如果超过自己的楼层在上升的楼层,每一层都会产生一个怒气值,比如,想在5层楼停,电梯却在8楼停那么6,7两层都会产生怒气值,然后在8楼下电梯走到5楼。

状态转移方程:A(i)=minj<i A(j)+Pi−1 k=j(i−k)∗sk +Pn k=i sk (j is the previous stop)

#include<bits/stdc++.h>
using namespace std;
const int maxn=1500+10;
const int inf=1e9+7;
int anger[maxn],s[maxn];
int main()
{
    freopen("in.txt","r",stdin);
    int t;
    int n;
    scanf("%d",&t);
    while(t--)
    {
      int sum=0;
      scanf("%d",&n);
      for(int i=1;i<=n;i++)
        {
            scanf("%d",&s[i]);
            sum+=s[i];
        }
      anger[0]=0;//0层大家上电梯的那一层愤怒值为0;
      for(int i=1;i<=n;i++)
      {
          sum-=s[i];//去往更高层的同学数量
          int skipanger=0;
          anger[i]=inf;
          for(int j=i-1;j>=0;j--)
        {
              //int k=j;
              anger[i]=min(anger[i],anger[j]+sum+skipanger);
              skipanger+=(i-j)*s[j];//去往底层学生的愤怒值
        }
      }
      cout<<anger[n]<<endl;
    }
    return 0;
}

clc;clear; % Driver script for solving the 1D Euler equations % Polynomial order used for approximation N = 6; Np = N + 1; % 每元素节点数 Nfaces = 2; % 1D元素面数 Nfp = 1; % 每个面的节点数 NODETOL = 1e-10; % 定义网格范围 xmin = 0.0; xmax = 1.0; K = 250; % 元素数量 % Generate simple mesh [Nv, VX, K, EToV] = MeshGen1D(xmin, xmax, K); % 构建参考元素 r = JacobiGL(0, 0, N); V = Vandermonde1D(N, r); invV = inv(V); Dr = Dmatrix1D(N, r, V); % 创建表面积分项 LIFT = Lift1D(Np, Nfaces, Nfp, V); % 构建所有节点的坐标 va = EToV(:,1)'; vb = EToV(:,2)'; x = ones(Np,1)*VX(va) + 0.5*(r+1)*(VX(vb)-VX(va)); % 计算几何因子 [rx, J] = GeometricFactors1D(x, Dr); % 计算边缘节点掩码 fmask1 = find( abs(r+1) < NODETOL)'; fmask2 = find( abs(r-1) < NODETOL)'; Fmask = [fmask1; fmask2]'; Fx = x(Fmask(:), :); % 构建表面法线和表面处的逆度量 [nx] = Normals1D(Nfp, Nfaces, K); % 构建连接矩阵 [EToE, EToF] = Connect1D(EToV); % 构建连接映射 [vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO] = BuildMaps1D(EToE, EToF, x, Np, Nfp, Nfaces, NODETOL, K, Fmask); % 修改返回 % 计算 Fscale Fscale = 1./(J(Fmask,:)); % 设置初始条件 -- Sod's problem gamma = 1.4; MassMatrix = inv(V')/V; cx = ones(Np,1)*sum(MassMatrix*x,1)/2; rho = ones(Np,K).*( (cx<0.5) + 0.125*(cx>=0.5)); rhou = zeros(Np,K); Ener = ones(Np,K).*((cx<0.5) + 0.1*(cx>=0.5))/(gamma-1.0); FinalTime = 0.2; % Solve Problem [rho,rhou,Ener] = Euler1D(rho, rhou, Ener, FinalTime, N, Np, K, r, x, VX, Dr, LIFT, nx, Fx, Fscale, vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO, Fmask, rx, J, gamma, V, invV, Nfp, Nfaces); % 添 x 参数 figure; plot(x, Ener); drawnow; pause(0.1); figure; subplot(3,1,1); plot(x, rho); title('Density'); subplot(3,1,2); plot(x, rhou./rho); title('Velocity'); subplot(3,1,3); plot(x, (gamma-1)*(Ener-0.5*rhou.^2./rho)); title('Pressure'); function [Nv, VX, K, EToV] = MeshGen1D(xmin, xmax, K) Nv = K+1; % Generate node coordinates VX = (1:Nv); for i = 1:Nv VX(i) = (xmax-xmin)*(i-1)/(Nv-1) + xmin; end % read element to node connectivity EToV = zeros(K, 2); for k = 1:K EToV(k,1) = k; EToV(k,2) = k+1; end end function [x] = JacobiGL(alpha, beta, N) x = zeros(N+1,1); if (N==1) x(1)=-1.0; x(2)=1.0; return; end [xint,w] = JacobiGQ(alpha+1, beta+1, N-2); x = [-1, xint', 1]'; end function [V1D] = Vandermonde1D(N, r) V1D = zeros(length(r), N+1); for j=1:N+1 V1D(:,j) = JacobiP(r(:), 0, 0, j-1); end end function [Dr] = Dmatrix1D(N, r, V) Vr = GradVandermonde1D(N, r); Dr = Vr/V; end function [LIFT] = Lift1D(Np, Nfaces, Nfp, V) Emat = zeros(Np, Nfaces*Nfp); % Define Emat Emat(1,1) = 1.0; Emat(Np,2) = 1.0; % inv(mass matrix)*\s_n (L_i,L_j)_{edge_n} LIFT = V*(V'*Emat); end function [rx, J] = GeometricFactors1D(x, Dr) xr = Dr*x; J = xr; rx = 1./J; end function [nx] = Normals1D(Nfp, Nfaces, K) nx = zeros(Nfp*Nfaces, K); % Define outward normals nx(1, :) = -1.0; nx(2, :) = 1.0; end function [EToE, EToF] = Connect1D(EToV) Nfaces = 2; % Find number of elements and vertices K = size(EToV,1); TotalFaces = Nfaces*K; Nv = K+1; % List of local face to local vertex connections vn = [1,2]; % Build global face to node sparse array SpFToV = spalloc(TotalFaces, Nv, 2*TotalFaces); sk = 1; for k=1:K for face=1:Nfaces SpFToV( sk, EToV(k, vn(face))) = 1; sk = sk+1; end end % Build global face to global face sparse array SpFToF = SpFToV*SpFToV' - speye(TotalFaces); % Find complete face to face connections [faces1, faces2] = find(SpFToF==1); % Convert face global number to element and face numbers element1 = floor( (faces1-1)/Nfaces ) + 1; face1 = mod( (faces1-1), Nfaces ) + 1; element2 = floor( (faces2-1)/Nfaces ) + 1; face2 = mod( (faces2-1), Nfaces ) + 1; % Rearrange into Nelements x Nfaces sized arrays ind = sub2ind([K, Nfaces], element1, face1); EToE = (1:K)'*ones(1,Nfaces); EToF = ones(K,1)*(1:Nfaces); EToE(ind) = element2; EToF(ind) = face2; end function [vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO] = BuildMaps1D(EToE, EToF, x, Np, Nfp, Nfaces, NODETOL, K, Fmask) % 修改函数定义 % number volume nodes consecutively nodeids = reshape(1:K*Np, Np, K); vmapM = zeros(Nfp, Nfaces, K); vmapP = zeros(Nfp, Nfaces, K); for k1=1:K for f1=1:Nfaces % find index of face nodes with respect to volume node ordering vmapM(:,f1,k1) = nodeids(Fmask(:,f1), k1); end end for k1=1:K for f1=1:Nfaces % find neighbor k2 = EToE(k1,f1); f2 = EToF(k1,f1); % find volume node numbers of left and right nodes vidM = vmapM(:,f1,k1); vidP = vmapM(:,f2,k2); x1 = x(vidM); x2 = x(vidP); % Compute distance matrix D = (x1 -x2 ).^2; if (D<NODETOL) vmapP(:,f1,k1) = vidP; end end end vmapP = vmapP(:); vmapM = vmapM(:); % Create list of boundary nodes mapB = find(vmapP==vmapM); vmapB = vmapM(mapB); % Create specific left (inflow) and right (outflow) maps mapI = 1; mapO = K*Nfaces; vmapI = 1; vmapO = K*Np; % 定义 vmapI, mapI, vmapO, mapO end function [rho, rhou, Ener] = Euler1D(rho, rhou, Ener, FinalTime, N, Np, K, r, x, VX, Dr, LIFT, nx, Fx, Fscale, vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO, Fmask, rx, J, gamma, V, invV, Nfp, Nfaces) % 添 x 参数 % Parameters CFL = 1.0; time = 0; % Prepare for adaptive time stepping mindx = min(x(2,:)-x(1,:)); % Limit initial solution rho = SlopeLimitN(rho, Np, V, invV, K, x, Dr); % 添 x 和 Dr 参数 rhou = SlopeLimitN(rhou, Np, V, invV, K, x, Dr); % 添 x 和 Dr 参数 Ener = SlopeLimitN(Ener, Np, V, invV, K, x, Dr); % 添 x 和 Dr 参数 % outer time step loop while(time < FinalTime) Temp = (Ener - 0.5*(rhou).^2./rho)./rho; cvel = sqrt(gamma*(gamma-1)*Temp); dt = CFL*min(min(mindx./(abs(rhou./rho)+cvel))); if(time+dt > FinalTime) dt = FinalTime-time; end % 3rd order SSP Runge-Kutta % SSP RK Stage 1. [rhsrho, rhsrhou, rhsEner] = EulerRHS1D(rho, rhou, Ener, N, Np, K, r, x, VX, Dr, LIFT, nx, Fx, Fscale, vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO, Fmask, rx, J, gamma, Nfp, Nfaces); % 添 vmapI, mapI, vmapO, mapO rho1 = rho + dt*rhsrho; rhou1 = rhou + dt*rhsrhou; Ener1 = Ener + dt*rhsEner; % Limit fields rho1 = SlopeLimitN(rho1, Np, V, invV, K, x, Dr); % 添 x 和 Dr 参数 rhou1 = SlopeLimitN(rhou1, Np, V, invV, K, x, Dr); % 添 x 和 Dr 参数 Ener1 = SlopeLimitN(Ener1, Np, V, invV, K, x, Dr); % 添 x 和 Dr 参数 % SSP RK Stage 2. [rhsrho, rhsrhou, rhsEner] = EulerRHS1D(rho1, rhou1, Ener1, N, Np, K, r, x, VX, Dr, LIFT, nx, Fx, Fscale, vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO, Fmask, rx, J, gamma, Nfp, Nfaces); % 添 vmapI, mapI, vmapO, mapO rho2 = (3*rho + rho1 + dt*rhsrho )/4; rhou2 = (3*rhou + rhou1 + dt*rhsrhou)/4; Ener2 = (3*Ener + Ener1 + dt*rhsEner)/4; % Limit fields rho2 = SlopeLimitN(rho2, Np, V, invV, K, x, Dr); % 添 x 和 Dr 参数 rhou2 = SlopeLimitN(rhou2, Np, V, invV, K, x, Dr); % 添 x 和 Dr 参数 Ener2 = SlopeLimitN(Ener2, Np, V, invV, K, x, Dr); % 添 x 和 Dr 参数 % SSP RK Stage 3. [rhsrho, rhsrhou, rhsEner] = EulerRHS1D(rho2, rhou2, Ener2, N, Np, K, r, x, VX, Dr, LIFT, nx, Fx, Fscale, vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO, Fmask, rx, J, gamma, Nfp, Nfaces); % 添 vmapI, mapI, vmapO, mapO rho = (rho + 2*rho2 + 2*dt*rhsrho )/3; rhou = (rhou + 2*rhou2 + 2*dt*rhsrhou)/3; Ener = (Ener + 2*Ener2 + 2*dt*rhsEner)/3; % Limit solution rho = SlopeLimitN(rho, Np, V, invV, K, x, Dr); % 添 x 和 Dr 参数 rhou = SlopeLimitN(rhou, Np, V, invV, K, x, Dr); % 添 x 和 Dr 参数 Ener = SlopeLimitN(Ener, Np, V, invV, K, x, Dr); % 添 x 和 Dr 参数 % Increment time and adapt timestep time = time + dt; end end function [rhsrho, rhsrhou, rhsEner] = EulerRHS1D(rho, rhou, Ener, N, Np, K, r, x, VX, Dr, LIFT, nx, Fx, Fscale, vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO, Fmask, rx, J, gamma, Nfp, Nfaces) % 修改函数定义 % compute maximum velocity for LF flux pres = (gamma-1.0)*(Ener - 0.5*(rhou).^2./rho); cvel = sqrt(gamma*pres./rho); lm = abs(rhou./rho)+cvel; % Compute fluxes rhof = rhou; rhouf = rhou.^2./rho + pres; Enerf = (Ener + pres).*rhou./rho; % Compute jumps at internal faces drho = zeros(Nfp*Nfaces, K); drho(:) = rho(vmapM) - rho(vmapP); drhou = zeros(Nfp*Nfaces, K); drhou(:) = rhou(vmapM) - rhou(vmapP); dEner = zeros(Nfp*Nfaces, K); dEner(:) = Ener(vmapM) - Ener(vmapP); drhof = zeros(Nfp*Nfaces, K); drhof(:) = rhof(vmapM) - rhof(vmapP); drhouf = zeros(Nfp*Nfaces, K); drhouf(:) = rhouf(vmapM) - rhouf(vmapP); dEnerf = zeros(Nfp*Nfaces, K); dEnerf(:) = Enerf(vmapM) - Enerf(vmapP); LFc = zeros(Nfp*Nfaces, K); LFc(:) = max(lm(vmapP), lm(vmapM)); % Compute fluxes at interfaces drhof(:) = nx(:).*drhof(:)/2.0 - LFc(:)/2.0.*drho(:); drhouf(:) = nx(:).*drhouf(:)/2.0 - LFc(:)/2.0.*drhou(:); dEnerf(:) = nx(:).*dEnerf(:)/2.0 - LFc(:)/2.0.*dEner(:); % Boundary conditions for Sod's problem rhoin = 1.000; rhouin = 0.0; pin = 1.000; Enerin = pin/(gamma-1.0); rhoout = 0.125; rhouout = 0.0; pout = 0.100; Enerout = pout/(gamma-1.0); % Set fluxes at inflow/outflow rhofin = rhouin; rhoufin = rhouin.^2./rhoin + pin; Enerfin = (pin/(gamma-1.0) + 0.5*rhouin^2/rhoin + pin).*rhouin./rhoin; lmI = lm(vmapI)/2; nxI = nx(mapI); % 使用 vmapI 和 mapI drhof (mapI) = nxI*(rhof (vmapI) - rhofin )/2.0 - lmI*(rho(vmapI) - rhoin); drhouf(mapI) = nxI*(rhouf(vmapI) - rhoufin)/2.0 - lmI*(rhou(vmapI) - rhouin); dEnerf(mapI) = nxI*(Enerf(vmapI) - Enerfin)/2.0 - lmI*(Ener(vmapI) - Enerin); rhofout = rhouout; rhoufout = rhouout.^2./rhoout + pout; Enerfout = (pout/(gamma-1.0) + 0.5*rhouout^2/rhoout + pout).*rhouout./rhoout; lmO = lm(vmapO)/2; nxO = nx(mapO); % 使用 vmapO 和 mapO drhof (mapO) = nxO*(rhof(vmapO) - rhofout)/2.0 - lmO*(rho (vmapO) - rhoout); drhouf(mapO) = nxO*(rhouf(vmapO) - rhoufout)/2.0 - lmO*(rhou(vmapO) - rhouout); dEnerf(mapO) = nxO*(Enerf(vmapO) - Enerfout)/2.0 - lmO*(Ener(vmapO) - Enerout); % compute right hand sides of the PDE's rhsrho = -rx.*(Dr*rhof) + LIFT*(Fscale.*drhof); rhsrhou = -rx.*(Dr*rhouf) + LIFT*(Fscale.*drhouf); rhsEner = -rx.*(Dr*Enerf) + LIFT*(Fscale.*dEnerf); end function ulimit = SlopeLimitN(u, Np, V, invV, K, x, Dr) % 添 x 和 Dr 参数 % Compute cell averages uh = invV*u; uh(2:Np,:)=0; uavg = V*uh; v = uavg(1,:); % Apply slope limiter as needed. ulimit = u; eps0 = 1.0e-8; % find end values of each element ue1 = u(1,:); ue2 = u(end,:); % find cell averages vk = v; vkm1 = [v(1), v(1:K-1)]; vkp1 = [v(2:K), v(K)]; % Apply reconstruction to find elements in need of limiting ve1 = vk - minmod([(vk-ue1); vk-vkm1; vkp1-vk]); ve2 = vk + minmod([(ue2-vk); vk-vkm1; vkp1-vk]); ids = find(abs(ve1-ue1) > eps0 | abs(ve2-ue2) > eps0); % Check to see if any elements require limiting if(~isempty(ids)) % create piecewise linear solution for limiting on specified elements uhl = invV*u(:,ids); uhl(3:Np,:)=0; ul = V*uhl; % apply slope limiter to selected elements ulimit(:,ids) = SlopeLimitLin(ul, x(:,ids), vkm1(ids), vk(ids), vkp1(ids), Np, Dr); % 添 Dr 参数 end end function ulimit = SlopeLimitLin(ul, xl, vm1, v0, vp1, Np, Dr) % 添 Dr 参数 % Compute various geometric measures ulimit = ul; h = xl(Np,:) - xl(1,:); x0 = ones(Np,1)*(xl(1,:) + h/2); hN = ones(Np,1)*h; % Limit function ux = (2./hN).*(Dr*ul); % 使用传入的 Dr ulimit = ones(Np,1)*v0 + (xl-x0).*(ones(Np,1)*minmod([ux(1,:); (vp1-v0)./h; (v0-vm1)./h])); end function mfunc = minmod(v) m = size(v,1); mfunc = zeros(1, size(v,2)); s = sum(sign(v),1)/m; ids = find(abs(s) == 1); if(~isempty(ids)) mfunc(ids) = s(ids).*min(abs(v(:,ids)), [], 1); end end function [x, w] = JacobiGQ(alpha, beta, N) if (N == 0) x(1) = -(alpha-beta)/(alpha+beta+2); w(1) = 2; return; end % Form symmetric matrix from recurrence. J = zeros(N+1); h1 = 2*(0:N) + alpha + beta; J = diag(-1/2*(alpha^2-beta^2)./(h1+2)./h1) + ... diag(2./(h1(1:N)+2).*sqrt((1:N).*((1:N)+alpha+beta).*... ((1:N)+alpha).*((1:N)+beta)./(h1(1:N)+1)./(h1(1:N)+3)),1); if (alpha+beta < 10*eps) J(1,1) = 0.0; end J = J + J'; % Compute quadrature by eigenvalue solve [V, D] = eig(J); x = diag(D); w = (V(1,:)').^2 * 2^(alpha+beta+1)/(alpha+beta+1) * gamma(alpha+1) * ... gamma(beta+1)/gamma(alpha+beta+1); end function [P] = JacobiP(x, alpha, beta, N) xp = x; dims = size(xp); if (dims(2) == 1) xp = xp'; end PL = zeros(N+1, length(xp)); % Initial values P_0(x) and P_1(x) gamma0 = 2^(alpha+beta+1)/(alpha+beta+1) * gamma(alpha+1) * ... gamma(beta+1)/gamma(alpha+beta+1); PL(1,:) = 1.0/sqrt(gamma0); if (N == 0) P = PL'; return; end gamma1 = (alpha+1)*(beta+1)/(alpha+beta+3)*gamma0; PL(2,:) = ((alpha+beta+2)*xp/2 + (alpha-beta)/2)/sqrt(gamma1); if (N == 1) P = PL(N+1,:)'; return; end % Repeat value in recurrence. aold = 2/(2+alpha+beta)*sqrt((alpha+1)*(beta+1)/(alpha+beta+3)); % Forward recurrence using the symmetry of the recurrence. for i=1:N-1 h1 = 2*i + alpha + beta; anew = 2/(h1+2)*sqrt( (i+1)*(i+1+alpha+beta)*(i+1+alpha)*... (i+1+beta)/(h1+1)/(h1+3)); bnew = - (alpha^2-beta^2)/h1/(h1+2); PL(i+2,:) = 1/anew*( -aold*PL(i,:) + (xp-bnew).*PL(i+1,:)); aold = anew; end P = PL(N+1,:)'; end function [DVr] = GradVandermonde1D(N, r) DVr = zeros(length(r), (N+1)); % Initialize matrix for i=0:N [DVr(:,i+1)] = GradJacobiP(r(:), 0, 0, i); end end function [dP] = GradJacobiP(r, alpha, beta, N) dP = zeros(length(r), 1); if(N == 0) dP(:,:) = 0.0; else dP = sqrt(N*(N+alpha+beta+1))*JacobiP(r(:), alpha+1, beta+1, N-1); end end 对上述代码分析,给出相对应的数学理论公式
07-10
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值